#ifndef _CENTERLINE_PANEL_H_
#define _CENTERLINE_PANEL_H_

#include "CenterLineThread.h"

#include "Common.h"
#include "Source.h"
#include "VolumeSource.h"
#include "ui_CenterLinePanel.h"
#include "ControlPanelBase.h"
#include "ControlPanelDock.h"
#include "DisplayModel.h"
#include "PipelineModel.h"
#include "SurfaceSource.h"
#include "DisplayFramework.h"
#include "View3D.h"
#include "FilteringThreadBase.h"
#include "SegmentationThreads.h"

#include "itkRegionOfInterestImageFilter.h"
#include "itkImageMaskSpatialObject.h"
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkSignedDanielssonDistanceMapImageFilter.h"
#include "itkApproximateSignedDistanceMapImageFilter.h"
#include "itkSignedMaurerDistanceMapImageFilter.h"
#include "itkAntiAliasBinaryImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkCastImageFilter.h"

#include <QDebug>
#include <QVector>
#include <QVector3D>
#include <QVector4D>
#include <vector>
#include <QStack>
#include <QQueue>

#include "vtkPropPicker.h"
#include "vtkSmartPointer.h"
#include "vtkSelectPolyData.h"
#include "vtkActor.h"
#include "vtkPoints.h"
#include "vtkLODActor.h"
#include "vtkClipPolyData.h"

//clips #include "clipscpp.h"


class PickingCallBackForVesselSeparation : public vtkCommand
{
protected:
	PickingCallBackForVesselSeparation()
	{
		m_pickedActor = NULL;
		m_displayModel = NULL;

		spacing[0] = 1; spacing[1]=1; spacing[2]=1;
		combox = NULL;
		keyControlPressed = false;
	}

public:
	static PickingCallBackForVesselSeparation * New()
	{
		return new PickingCallBackForVesselSeparation;
	}

	void Delete()
	{
		delete this;
	}

	virtual void Execute(vtkObject *caller, unsigned long eventId, void *callData)
	{
		vtkPropPicker* picker = reinterpret_cast<vtkPropPicker*>(caller);
		vtkActor* currentActor = picker->GetActor();

		// picked nothing
		if ( currentActor == NULL)
		{
			cout << "picked nothing" << endl;
			return;
		}


		// previous actor
		if ( m_pickedActor)
		{
			if ( m_displayModel->isASphere(m_pickedActor))
				m_pickedActor->GetProperty()->SetColor(previousColor);	// set to white
			else if ( m_displayModel->isABranch(m_pickedActor))
			{
				if ( keyControlPressed == false)
					m_displayModel->clearSelectedBranches();
				else
					;
					//m_displayModel->setSelectABranch(m_pickedActor, false);
			}
		}

		if ( m_pickedActor == picker->GetActor())	//deselect
		{
			cout << "picked the same actor" << endl;
			m_pickedActor = NULL;
			return;
		}

		double *pos = picker->GetPickPosition();
		int transfered[4] = {0};
		transfered[0] = int(pos[0]/spacing[0]);
		transfered[1] = int(pos[1]/spacing[1]);
		transfered[2] = int(pos[2]/spacing[2]);

		m_pickedActor = picker->GetActor();

		if ( m_displayModel->isASphere(m_pickedActor))
		{
			double* icolor = m_pickedActor->GetProperty()->GetColor();
			previousColor[0] = icolor[0];
			previousColor[1] = icolor[1];
			previousColor[2] = icolor[2];

			if ( combox != NULL)
			{
				for (int i = 0; i < combox->count(); i++)
				{
					QStringList sl = (combox->itemText(i)).split(" ");
					if (std::labs(sl[0].toUInt()-transfered[0]) <= 5 && 
						std::labs(sl[1].toUInt()-transfered[1])<=5 && 
						std::labs(sl[2].toUInt()-transfered[2])<=5)
					{
						combox->setCurrentIndex(i);
						cout << "transfered:" << transfered[0] << ", " << transfered[1] << ", " << transfered[2] << ", " << sl[3].toDouble()<< endl;
						cout << "matched:" << sl[0].toInt() << ", " << sl[1].toInt() << ", " << sl[2].toInt() << ", " << sl[3].toDouble() << endl;
						break;
					}
				}
			}

			m_pickedActor->GetProperty()->SetColor(1, 0, 0);
		}
		else if ( m_displayModel->isALine(m_pickedActor))
		{
			m_pickedActor->GetProperty()->SetColor(1, 0, 0);
		}
		else if ( m_displayModel->isABranch(m_pickedActor))
		{
			cout << "picked a branch" << endl;

			m_displayModel->setSelectABranch(m_pickedActor, true);
		}
		else
		{
			cout << "picked unknown actor" << endl;
			return;
		}

		m_displayModel->updateActive3DScene();
	}

	QString getPickedPosition()
	{
		if ( m_pickedActor)
			return m_displayModel->getSpherePosition(m_pickedActor);

		return QString("");
	}

	vtkActor* getPickedActor()
	{
		return m_pickedActor;
	}

public:
	DisplayModel* m_displayModel;
	double		spacing[3];
	QComboBox	*combox;

	//QList<

	QList<vtkSmartPointer<vtkActor>> m_pickedBranchActors;
	vtkSmartPointer<vtkActor> m_pickedActor;
	double	previousColor[3];
	bool	keyControlPressed;
};

typedef struct _branchProperty
{
	bool extreme;				//it's extreme point
	bool useless;				//this branch will be deleted
}BranchProperty_t;

#define NODE_TYPE_PORTAL	0
#define NODE_TYPE_HEPATIC	1


typedef struct _branchPoint
{
	int id;
	_branchPoint* Parent;
	QList<_branchPoint*> Children;
	bool readyToRemove;
	bool removed;
	int type;
	int maxDepth;
	int depthFromPortal;
	int subNodeNum;

	int group;

	QVector4D RootPoint;
	QList<QList<QVector4D>*> SubBranches;

	QList<BranchProperty_t*> SubBranchProperties;
}BranchNode_t;


class CenterLinePanel : public ControlPanelBase
{
	Q_OBJECT

	typedef itk::Image<float, 3> CImageTypeFloat3D;

	PickingCallBackForVesselSeparation *m_pickingCallBack;


private:
	Ui::CenterLinePanel m_ui;
	QString m_currentMaskFile;

	Source *m_sourceSeed;
	Source *m_sourceSkeleton;
	Source *m_sourceLiver;
	double* m_currentSpacing;

	View3D *m_ballDisplayingView;
	QMap<QString, QVector<int> > m_seedPoints;

	QStringList m_portalPoint;
	QList<BranchNode_t*> m_branchNodeList;

	int		m_iSizeSphere; 
	int*	m_pIndexSphereX;
	int*	m_pIndexSphereY;
	int*	m_pIndexSphereZ;


public:
	CenterLinePanel(QWidget *parent);

	void setSource(const QModelIndex& index);
	void setDisplayModel(DisplayModel* dm);
	void removeCurrentSelectedActor();

private:
	QString transferUI2Display(const QString& text);	//transfer 
	void processDivisionByClick();
	void divideLiverByClicking();
	SurfaceSource* volume2SmoothSurface(Source *source, QString label, int intenal);

private slots:
	void slotExtractCenterLine();
	void slotLoadMask();
	void slotProcessMask();
	void slotLiverDivision2();
	void slotLiverDivision4(){}
	void slotPointsSelected();
	void slotSelectSeed();
	void slotCenterLineProgressDone();

	void slotSaveBranchPoints();
	void slotLoadBranchPoints();
	void slotShowAllPoints();
	void slotClearAllPoints();
	void slotLocateHepaticVein();
	void slotSliderUpperMoved(int);
	void slotSliderLowerMoved(int);

	void slotProcessKeyEventFromView(QKeyEvent* event);

	void slotRebuildVessel();
	void slotRemoveVessel();
	void slotRemoveHepatic();

	void slotLocatePortalPoint();
	void slotAutoClassifyBranches();
	void slotAutoLeftLiver();
	void slotAutoRightLiver();
	void slotSetSkeletonSource();
	void slotSetLiverSource();
	void slotSmoothVolume();
	void slotCreateSmoothSurface();

	void slotTesting1();
	void slotTesting2();
	void slotTesting3();
	void slotTestClips();


private:
	FilteringThreadBase *m_threadCenterLineExtraction;
	std::vector<Node> m_vectorControlPoints;

	void addSphere(const QString &text);
	void removeSphere(const QString &text);
	double transferRadius(double r);
	CImageTypeFloat3D::Pointer m_pDistanceMap;


	template<class TPixelType>
	void createCenterLineThread()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		m_threadCenterLineExtraction = new CenterLineThread<TPixelType>(currentSource());
		((CenterLineThread<TPixelType>*)m_threadCenterLineExtraction)->SetBiggestBranch(m_ui.m_cbBiggestBranch->isChecked());

		connect(m_threadCenterLineExtraction, SIGNAL(emitProgressStart(QString)), m_parent, SLOT(slotProgressStart(QString)));
		connect(m_threadCenterLineExtraction, SIGNAL(emitProgressUpdate(int)), m_parent, SLOT(slotProgressUpdate(int)));
		connect(m_threadCenterLineExtraction, SIGNAL(emitProgressDone()), m_parent, SLOT(slotProgressDone()));
		connect(m_threadCenterLineExtraction, SIGNAL(emitProgressDone()), SLOT(slotCenterLineProgressDone()));

		m_threadCenterLineExtraction->start();
	}


	template<class TPixelType>
	double* getCurrentSpacing()
	{
		typedef itk::Image<TPixelType, 3> TImageType;

		double *spacing = new double[3];

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename TImageType::SpacingType spacingtype = m_inputImage->GetSpacing();

		spacing[0] = spacingtype[0];
		spacing[1] = spacingtype[1];
		spacing[2] = spacingtype[2];

		return spacing;
	}

	template<class TPixelType>
	void createMaskedVolume()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::ImageFileReader<TImageType> ImageFileReaderType;

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename ImageFileReaderType::Pointer reader = ImageFileReaderType::New();
		typename TImageType::Pointer pVTKMaskFile;

		typename ImageDuplicator::Pointer dup = ImageDuplicator::New();
		dup->SetInputImage(m_inputImage);
		dup->Update();
		dup->GetOutput()->FillBuffer(0);
		typename TImageType::Pointer m_dupImage = dup->GetOutput();

		typename TImageType::SizeType size = m_inputImage->GetBufferedRegion().GetSize();
		int bytes = size[0]*size[1]*size[2];

		//RxMaskVolume *pMaskVolume = new RxMaskVolume();
		//pMaskVolume->CreateMaskVolume(size[0], size[1], size[2]);

		FILE* fp;
		unsigned char * maskbuf = new unsigned char[bytes];
		bool isVTKFile = false;
		if ( m_currentMaskFile.at(m_currentMaskFile.count()-1) == QChar('k'))	//vtk file
		{
			reader->SetFileName(m_currentMaskFile.toLatin1().data());
			reader->Update();
			delete [] maskbuf;
			isVTKFile = true;
			pVTKMaskFile = reader->GetOutput();
		}
		else
		{
			fopen_s(&fp, m_currentMaskFile.toLatin1().data(), "rb");
			//fread(pMaskVolume->m_pbyMask, sizeof(unsigned char), pMaskVolume->m_iMaskSize*8, fp);
			fread(maskbuf, sizeof(unsigned char), bytes, fp);
			fclose(fp);
		}

		typename TImageType::IndexType pixelIndex;
		for (unsigned int i = 0; i < size[2]; i++)
		{
			for (unsigned int j = 0; j< size[1]; j++ )
			{
				for (unsigned int k = 0; k < size[0]; k++)
				{
					pixelIndex[0] = k;
					pixelIndex[1] = j;
					pixelIndex[2] = i;//+m_ui.m_spinMaskOffset->value();
					if ( pixelIndex[2] >= size[2])
						pixelIndex[2] = i;

					int location = i*size[0]*size[1] + j * size[0] + k;

					if ( isVTKFile)
					{
						if ( pVTKMaskFile->GetPixel(pixelIndex) == 0)
						{
							m_inputImage->SetPixel(pixelIndex, 0);
						}
						else
							m_dupImage->SetPixel(pixelIndex, 1);
					}
					else
					{
						if ( maskbuf[location] == 0)
						{
							m_inputImage->SetPixel(pixelIndex, 0);
						}
						else
							m_dupImage->SetPixel(pixelIndex, 1);
					}



					/*if ( pMaskVolume->GetMask(k,j,i) == 0)
						m_inputImage->SetPixel(pixelIndex, 0);
					else
						m_dupImage->SetPixel(pixelIndex, 1);*/
				}
			}
		}


		VolumeSource<TImageType> *subsource2 = 
			new VolumeSource<TImageType>(QString("MaskVolume#%1").arg(currentSource()->rowCount()+1));
		subsource2->setVolume(m_dupImage);
		currentSource()->addChild(subsource2);

		VolumeSource<TImageType> *subsource = 
			new VolumeSource<TImageType>(QString("MaskedVolume#%1").arg(currentSource()->rowCount()+1));
		subsource->setVolume(m_inputImage);
		currentSource()->addChild(subsource);
	}



	template<class TPixelType>
	void processSelectedPoints()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();

		typename ImageDuplicator::Pointer dup = ImageDuplicator::New();
		dup->SetInputImage(m_inputImage);
		dup->Update();
		dup->GetOutput()->FillBuffer(0);

		typename TImageType::Pointer m_dupVolume = dup->GetOutput();

		typename TImageType::IndexType pixelIndex;
		typename TImageType::IndexType averageIndex;
		typename TImageType::SizeType size;
		typename TImageType::SpacingType spacing;

		averageIndex[0] = 0; averageIndex[1] = 0; averageIndex[2] = 0;
		size = m_inputImage->GetBufferedRegion().GetSize();
		spacing = m_inputImage->GetSpacing();

		vtkDataSet *dataSet  = m_displayModel->getSelectedPoints();
		int numberOfPoints = dataSet->GetNumberOfPoints();
		for (vtkIdType i = 0; i < numberOfPoints; i++)
		{
			double *data = dataSet->GetPoint(i);

			pixelIndex[0] = data[0] / spacing[0];
			pixelIndex[1] = data[1] / spacing[1];
			pixelIndex[2] = data[2] / spacing[2];

			//std::cout << "[" << i << "]:" << pixelIndex[0] << ", " << pixelIndex[1] << ", " << pixelIndex[2] << std::endl;
			m_dupVolume->SetPixel(pixelIndex, 1);
		}

		double* bounds = dataSet->GetBounds();

		averageIndex[0] = ((bounds[1] + bounds[0])/2)/spacing[0] + 1;
		averageIndex[1] = ((bounds[3] + bounds[2])/2)/spacing[1] + 1;
		averageIndex[2] = ((bounds[5] + bounds[4])/2)/spacing[2] + 1;


		static int pointNum = 1;
		QString str = QString("PointVolume#%1").arg(pointNum++);

		VolumeSource<TImageType> *subsource = new VolumeSource<TImageType>(str);
		QVector<int> vseed; 
		vseed.append((int)averageIndex[0]);
		vseed.append((int)averageIndex[1]);
		vseed.append((int)averageIndex[2]);
		vseed.append((int)(averageIndex[0]*spacing[0]));
		vseed.append((int)(averageIndex[1]*spacing[1]));
		vseed.append((int)(averageIndex[2]*spacing[2]));

		m_seedPoints.insert(str,vseed);
		subsource->setVolume(m_dupVolume);
		m_pipelineModel->addSource(subsource);

	}

	double Distance(QVector4D* v1, QVector4D* v2)
	{
		return sqrt(
			(v1->x()-v2->x())*(v1->x()-v2->x()) +
			(v1->y()-v2->y())*(v1->y()-v2->y()) +
			(v1->z()-v2->z())*(v1->z()-v2->z()));

	}
	void Normalize(QVector4D* v, double r)
	{
		double l;

		l = r / sqrt(v->x()*v->x() + v->y()*v->y() + v->z()*v->z());
		v->setX(v->x()*l);
		v->setY(v->y()*l);
		v->setZ(v->z()*l);
	}


	QVector4D* pickPointsOnSphere(int *center, double radius, const int pointNum)
	{
		QVector4D *vectors = new QVector4D[pointNum];

		//initial points
		for (int i = 0; i < pointNum; i++)
		{
			vectors[i].setX(rand()%1000 - 500);
			vectors[i].setY(rand()%1000 - 500);
			vectors[i].setZ(rand()%1000 - 500);

			Normalize(&vectors[i], radius);

			printf("random points %f,%f,%f\n", vectors[i].x(), vectors[i].y(), vectors[i].z());
		}


		// average the distribution

		int count = 0;
		int countMax = 10000;
		int minPointIndex1;		// index of array
		int minPointIndex2;		
		int maxPointIndex1;		// index of array
		int maxPointIndex2;		

		double minDist;
		double maxDist;
		double dist;
		QVector4D *point1;
		QVector4D *point2;
		while ( count++ < countMax)
		{
			minPointIndex1 = 0;
			minPointIndex2 = 1;
			maxPointIndex1 = 0;
			maxPointIndex2 = 1;
			maxDist = minDist = Distance(&vectors[minPointIndex1], &vectors[minPointIndex2]);

			for (int i = 0; i < pointNum - 1; i++)
			{
				for (int j = i+1; j < pointNum; j++)
				{
                    dist = Distance(&vectors[i], &vectors[j]);
					if ( dist < minDist)
					{
						minDist = dist; minPointIndex1 = i; minPointIndex2 = j;
					}
					
					if ( dist > maxDist)
					{
						maxDist = dist; maxPointIndex1 = i; maxPointIndex2 = j;
					}
				}
			}

			point1 = &vectors[minPointIndex1];
			point2 = &vectors[minPointIndex2];

			// enlarge the distance between two closest points
			vectors[minPointIndex2].setX(point1->x() + 1.2 * (point2->x()-point1->x()));
			vectors[minPointIndex2].setY(point1->y() + 1.2 * (point2->y()-point1->y()));
			vectors[minPointIndex2].setZ(point1->z() + 1.2 * (point2->z()-point1->z()));

			vectors[minPointIndex1].setX(point1->x() - 0.2 * (point2->x()-point1->x()));
			vectors[minPointIndex1].setY(point1->y() - 0.2 * (point2->y()-point1->y()));
			vectors[minPointIndex1].setZ(point1->z() - 0.2 * (point2->z()-point1->z()));

			Normalize(&vectors[minPointIndex1], radius);
			Normalize(&vectors[minPointIndex2], radius);


			// reduce the distance between two furthest points

// 			point1 = &vectors[maxPointIndex1];
// 			point2 = &vectors[maxPointIndex2];

// 			vectors[maxPointIndex1].setX(point1->x() + 0.1 * (point2->x()-point1->x()));
// 			vectors[maxPointIndex1].setY(point1->y() + 0.1 * (point2->y()-point1->y()));
// 			vectors[maxPointIndex1].setZ(point1->z() + 0.1 * (point2->z()-point1->z()));
// 
// 			vectors[maxPointIndex2].setX(point1->x() - 0.1 * (point2->x()-point1->x()));
// 			vectors[maxPointIndex2].setY(point1->y() - 0.1 * (point2->y()-point1->y()));
// 			vectors[maxPointIndex2].setZ(point1->z() - 0.1 * (point2->z()-point1->z()));
// 
// 			Normalize(&vectors[maxPointIndex1], radius);
// 			Normalize(&vectors[maxPointIndex2], radius);

		}

		return vectors;
	}
	
	template<class TPixelType>
	bool IsInVolume(int x, int y, int z)
	{
		typedef itk::Image<TPixelType, 3> TImageType;

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename TImageType::IndexType index;
		index[0] = x;
		index[1] = y;
		index[2] = z;

		return m_inputImage->GetPixel(index) == 0 ? false : true;
	}

	void FreeSphereStructure()
	{
		if (m_pIndexSphereX) {
			delete [] m_pIndexSphereX;
			m_pIndexSphereX = NULL;
		}
		if (m_pIndexSphereY) {
			delete [] m_pIndexSphereY;
			m_pIndexSphereY = NULL;
		}
		if (m_pIndexSphereZ) {
			delete [] m_pIndexSphereZ;
			m_pIndexSphereZ = NULL;
		}

		m_iSizeSphere = 0;
	}

	void MakeSphereStructure(int nMaxRadius)
	{
		int* pSphereBox = new int[(nMaxRadius + 1)*(nMaxRadius + 1)*(nMaxRadius + 1)];
		int x, y, z, iIndex = 0, indexCount = 0, iRadius, i, j, k;

		FreeSphereStructure();

		for (z = 0; z <= nMaxRadius; z++)
		{
			for (y = 0; y <= nMaxRadius; y++)
			{
				for (x = 0; x <= nMaxRadius; x++) 
				{
					pSphereBox[iIndex] = (int)(sqrt((double)(x*x + y*y + z*z)) + 0.5);
					if (0 < pSphereBox[iIndex] && pSphereBox[iIndex] <= nMaxRadius) 
					{
						if ((x != 0 && y == 0 && z == 0) || (x == 0 && y != 0 && z == 0) || (x == 0 && y == 0 && z != 0))									
							m_iSizeSphere += 2;
						if ((x == 0 && y != 0 && z != 0) || (x != 0 && y == 0 && z != 0) || (x != 0 && y != 0 && z == 0))									
							m_iSizeSphere += 4;
						if (x != 0 && y != 0 && z != 0)									
							m_iSizeSphere += 8;
					}
					iIndex++;
				}
			}
		}

		m_pIndexSphereX = new int[m_iSizeSphere];
		m_pIndexSphereY = new int[m_iSizeSphere];
		m_pIndexSphereZ = new int[m_iSizeSphere];

		memset(m_pIndexSphereX, 0, m_iSizeSphere*sizeof(int));
		memset(m_pIndexSphereY, 0, m_iSizeSphere*sizeof(int));
		memset(m_pIndexSphereZ, 0, m_iSizeSphere*sizeof(int));


		indexCount = 0;
		for (iRadius = nMaxRadius-1 ; iRadius <= nMaxRadius; iRadius++) 
		{
			iIndex = 0;
			for (z = 0; z <= nMaxRadius; z++)
			{
				for (y = 0; y <= nMaxRadius; y++)
				{
					for (x = 0; x <= nMaxRadius; x++)	
					{
						if (pSphereBox[iIndex] == iRadius) 
						{
							if (x != 0 && y == 0 && z == 0)
								for (i = -1; i <= 1; i+=2) 
								{
									m_pIndexSphereX[indexCount] = x*i;
									m_pIndexSphereY[indexCount] = 0;
									m_pIndexSphereZ[indexCount] = 0;
									indexCount++;
								}
							if (x == 0 && y != 0 && z == 0)
								for (j = -1; j <= 1; j+=2) 
								{
									m_pIndexSphereX[indexCount] = 0;
									m_pIndexSphereY[indexCount] = y*j;
									m_pIndexSphereZ[indexCount] = 0;
									indexCount++;
								}

							if (x == 0 && y == 0 && z != 0)
								for (k = -1; k <= 1; k+=2) 
								{
									m_pIndexSphereX[indexCount] = 0;
									m_pIndexSphereY[indexCount] = 0;
									m_pIndexSphereZ[indexCount] = z*k;
									indexCount++;
								}

							if (x == 0 && y != 0 && z != 0)
								for (k = -1; k <= 1; k+=2)
									for (j = -1; j <= 1; j+=2) 
									{
										m_pIndexSphereX[indexCount] = 0;
										m_pIndexSphereY[indexCount] = y*j;
										m_pIndexSphereZ[indexCount] = z*k;
										indexCount++;
									}

							if (x != 0 && y == 0 && z != 0)
								for (k = -1; k <= 1; k+=2)
									for (i = -1; i <= 1; i+=2) 
									{
										m_pIndexSphereX[indexCount] = x*i;
										m_pIndexSphereY[indexCount] = 0;
										m_pIndexSphereZ[indexCount] = z*k;
										indexCount++;
									}
							if (x != 0 && y != 0 && z == 0)
								for (j = -1; j <= 1; j+=2)
									for (i = -1; i <= 1; i+=2) 
									{
										m_pIndexSphereX[indexCount] = x*i;
										m_pIndexSphereY[indexCount] = y*j;
										m_pIndexSphereZ[indexCount] = 0;
										indexCount++;
									}


							if (x != 0 && y != 0 && z != 0) 
							{									
								for (k = -1; k <= 1; k+=2) 
								{
									for (j = -1; j <= 1; j+=2) 
									{
										for (i = -1; i <= 1; i+=2)
										{
											m_pIndexSphereX[indexCount] = x*i;
											m_pIndexSphereY[indexCount] = y*j;
											m_pIndexSphereZ[indexCount] = z*k;
											indexCount++;
										}
									}
								}
							}
						}
						iIndex++;
					}
				}
			}
		}

		delete [] pSphereBox;
	}




	template<class TPixelType>
	double evaluateBranchRadius(int *seed, double radius)
	{
		typedef itk::Image<TPixelType, 3> TImageType;

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename TImageType::Pointer m_parentImage = ((VolumeSource<TImageType>*)currentSource()->parent())->getVolume();
		if (!m_parentImage)
		{
			cout << "Distance map is missing, cannot perform separation" << endl;
			return -1;
		}

		typename TImageType::IndexType inputIndex;
		inputIndex[0] = seed[0];
		inputIndex[1] = seed[1];
		inputIndex[2] = seed[2];

		//printf("Evaluate Branch point [%d,%d,%d, %f]\n",seed[0],seed[1],seed[2],radius);

		int offset = 10;
		int maxExpand = 5;
        
		typename TImageType::IndexType finalIndex;
		typename TImageType::IndexType tempIndex;
		bool validRadius = true;
		double finalValue = radius;
		for (int i = seed[0]-offset; i < seed[0]+offset; i++)
		{
			for (int j = seed[1]-offset; j < seed[1]+offset; j++)
			{
				for ( int k = seed[2]-offset; k < seed[2]+offset; k++)
				{
					for (int r = (int)radius; r <= (int)radius+maxExpand; r+=1)
					{
						validRadius = true;
						MakeSphereStructure(r);
						int emptyCount = 0;
						for (int m = 0; m < m_iSizeSphere; m++)
						{
							tempIndex[0] = m_pIndexSphereX[m]+i;
							tempIndex[1] = m_pIndexSphereY[m]+j;
							tempIndex[2] = m_pIndexSphereZ[m]+k;

							//maybe too restrict
							if (m_parentImage->GetPixel(tempIndex) == 0)		
							{
								emptyCount++;
							}

							if ( emptyCount >= 10)
							{
								validRadius = false;
								break;
							}
						}

						//printf("Empty count [%d] in radius [%d]\n", emptyCount, r);

						if ( validRadius == true && r > finalValue)
						{
							finalValue = r;
							finalIndex[0] = i;
							finalIndex[1] = j;
							finalIndex[2] = k;
						}
					}
				}
			}
		}


		if ( finalValue != radius)
		{
			printf("Branch point [%d,%d,%d,%f] moved to [%d,%d,%d,%f]\n", (int)seed[0], (int)seed[1], (int)seed[2], (float)radius,
				(int)finalIndex[0], (int)finalIndex[1], (int)finalIndex[2], (float)finalValue);
		}
		else
		{
			printf("Branch point didn't change\n");
		}

		return finalValue;
	}


	template<class TPixelType>
	void removeABranch(BranchNode_t* branch)
	{
/*
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();

		typename ImageDuplicator::Pointer dup0 = ImageDuplicator::New();
		dup0->SetInputImage(m_inputImage);
		dup0->Update();
		typename TImageType::Pointer m_outputImage = dup0->GetOutput();

		TImageType::IndexType branchIndex;

		QList<QVector4D> removingPoints;
		QStack<BranchPoint_t*> branches;
		
		branches.push(branch);

		while(1)
		{
			if ( branches.count() == 0)
				break;

			BranchPoint_t *bp = branches.pop();

			for (int i = 0; i < bp->PointsList.count(); i++)
			{
				removingPoints.push_back(bp->PointsList.at(i));
			}

			for (int i = 0; i < bp->ChildBranchPoints.count(); i++)
			{
				for (int j = 0; j < branchPointList.count(); j++)
				{
					if (branchPointList.at(j)->RootPoint ==  bp->ChildBranchPoints.at(i))
					{
						branches.push(branchPointList.at(j));
					}
				}
			}
		}

		
		for (int i = 0; i < removingPoints.count(); i++)
		{
			branchIndex[0] = removingPoints.at(i).x();
			branchIndex[1] = removingPoints.at(i).y();
			branchIndex[2] = removingPoints.at(i).z();

			m_outputImage->SetPixel(branchIndex, 0);

			cout << "removed point [" << branchIndex[0] << ", " << branchIndex[1] << ", " << branchIndex[2] << endl;
		}


		VolumeSource<TImageType> *subsource = 
			new VolumeSource<TImageType>(QString("test#%1").arg(currentSource()->rowCount()+1));
		subsource->setVolume(m_outputImage);
		currentSource()->addChild(subsource);*/

	}


	void calculateMaxDepth()
	{
		int maxID = 0;
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			if ( m_branchNodeList[i]->id > maxID)
				maxID = m_branchNodeList[i]->id;
		}

		int* depth = new int[maxID+1];

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* bn = m_branchNodeList[i];
			memset(depth, 0, maxID*sizeof(int));

			QStack<BranchNode_t*> tmpStack;
			tmpStack.push(bn);

			while(1)
			{
				if ( tmpStack.isEmpty()) break;

				BranchNode_t* node = tmpStack.pop();

				for (int j = 0; j < node->Children.count(); j++)
				{
					tmpStack.push(node->Children[j]);
					depth[(node->Children[j])->id] = depth[node->id];
					(depth[(node->Children[j])->id])++;
				}
			}

			bn->maxDepth = 0;
			bn->subNodeNum = 0;
			for(int j = 0; j < maxID; j++)
			{
				if ( depth[j] > 0)
					bn->subNodeNum++;

				if ( depth[j]>bn->maxDepth)
					bn->maxDepth = depth[j];
			}

			qDebug() << bn->RootPoint << ", d=" << bn->maxDepth << "s=" << bn->subNodeNum;
		}
	}


	// this function make sure to choose a point in the portal vein
	QStringList autoLocatePortalPoint()
	{
		BranchNode_t *root = NULL;
		for (int i = 0; i < m_branchNodeList.size(); i++)
		{
			if ( m_branchNodeList[i]->Parent == NULL)
			{
				root = m_branchNodeList[i];
				break;
			}
		}
		
		if ( !root)
			return QStringList();

		QQueue<BranchNode_t*> tmpQueue;
		BranchNode_t* portalPoint = NULL;
		tmpQueue.enqueue(root);

		while (1)
		{
			if ( tmpQueue.isEmpty()) break;

			BranchNode_t *bn = tmpQueue.dequeue();

			if ( bn->Children.size() == 0)
				continue;

			int depthLargerThan4 = 0;
			QList<int> subNodeNums;
			for (int i = 0; i < bn->Children.size(); i++)
			{
				if (bn->Children[i]->maxDepth >= 4)
				{
					subNodeNums.push_back(bn->Children[i]->subNodeNum);
					depthLargerThan4++;
				}
				
				tmpQueue.enqueue(bn->Children[i]);
			}

			if ( depthLargerThan4 >= 2)
			{
				
				portalPoint = bn;
				break;
			}
		}

		if ( !portalPoint)
		{
			cout << "no portal point" << endl;
			return QStringList();
		}

		QString str = QString("%1 %2 %3 %4").arg(portalPoint->RootPoint.x())
			.arg(portalPoint->RootPoint.y())
			.arg(portalPoint->RootPoint.z())
			.arg(portalPoint->RootPoint.w());

		qDebug() << "got portal point " << str;

		return str.split(" ");
	}


	bool isValidNode(BranchNode_t* node, int childNum, int maxDepth)
	{
		int validChildren = 0;
		for (int i = 0; i < node->Children.size(); i++)
		{
			if (node->Children[i]->maxDepth >= maxDepth)
			{
				validChildren++;
			}
		}

		if ( validChildren >= childNum)
			return true;
		else
			return false;
	}


	void removeNode(BranchNode_t* node, int allChildren = true)
	{
		QQueue<BranchNode_t*> tmpQueue;
		BranchNode_t* targetNode = NULL;
		tmpQueue.enqueue(node);

		while (1)
		{
			if ( tmpQueue.isEmpty()) break;

			BranchNode_t *bn = tmpQueue.dequeue();

			for (int i = 0; i < bn->Children.count(); i++)
			{
				m_displayModel->removeABranch(bn->id, i);
				bn->SubBranchProperties[i]->useless = true;

				tmpQueue.enqueue(bn->Children[i]);
			}
		}
	}



	//find out a child, which has childnum children and for each has at least maxDepth
	BranchNode_t* getNextValidNode(BranchNode_t* node, int childNum, int maxDepth)
	{
		QQueue<BranchNode_t*> tmpQueue;
		BranchNode_t* targetNode = NULL;
		tmpQueue.enqueue(node);

		while (1)
		{
			if ( tmpQueue.isEmpty()) break;

			BranchNode_t *bn = tmpQueue.dequeue();

			if ( bn->Children.count() == 0)
				continue;

			qDebug() << "current node " << bn->RootPoint;

			int validChildren = 0;
			for (int i = 0; i < bn->Children.count(); i++)
			{
				if (bn->SubBranchProperties[i]->useless == true)
					continue;

				if (bn->Children[i]->maxDepth >= maxDepth) 
				{
					qDebug() << "valid--> " << bn->Children[i]->RootPoint;
					validChildren++;
				}

				tmpQueue.enqueue(bn->Children[i]);
			}

			if ( validChildren >= childNum)
			{
				targetNode = bn;
				break;
			}
		}

		return targetNode;
	}


	template<class TImageType>
	void selectGroupingBranches(QList<BranchNode_t*> &nodes)
	{
		int branch_color = 1;
		QList<int> group;

		for ( int i = 0; i < nodes.count(); i++)
		{
			if(nodes[i] == NULL)
				continue;

			QVector4D cvt_pnode = nodes[i]->Parent->RootPoint;
			cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
			cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
			cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

			QVector4D cvt_cnode =  nodes[i]->RootPoint;
			cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
			cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
			cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

			if ( m_ui.m_btnWithLiver->isChecked())
				m_displayModel->clearSelectedBranches();

			if ( nodes[i]->group == 0)
			{
				m_displayModel->setBranchClass(branch_color,  cvt_pnode, cvt_cnode);

				if ( m_ui.m_btnWithLiver->isChecked())
					divideBySelectedBranch<TImageType>(branch_color);
				branch_color++;
			}
			else
			{
				int exist = 0;
				for (int j = 0; j < group.count(); j++)
				{
					if ( group[j] == nodes[i]->group)
						exist = 1;
				}

				if ( exist == 0)
					group.push_back(nodes[i]->group);

			}

		}

		for (int i = 0; i < group.count(); i++)
		{
			if ( m_ui.m_btnWithLiver->isChecked())
				m_displayModel->clearSelectedBranches();

			for (int j = 0; j < nodes.count(); j++)
			{
				QVector4D cvt_pnode = nodes[j]->Parent->RootPoint;
				cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
				cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
				cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

				QVector4D cvt_cnode =  nodes[j]->RootPoint;
				cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
				cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
				cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

				if ( group[i] == nodes[j]->group)
					m_displayModel->setBranchClass(nodes[i]->group,  cvt_pnode, cvt_cnode);
			}

			if ( m_ui.m_btnWithLiver->isChecked())
				divideBySelectedBranch<TImageType>(nodes[i]->group);
		}


		//group.push_back(nodes[i]->group);
		
		//if ( m_ui.m_btnWithLiver->isChecked())		
			

	}


	template<class TPixelType>
	void autoClassifyBranches_bak_4parts()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		QStringList seedlist = autoLocatePortalPoint();
		if ( seedlist.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return;
		}


		typename TImageType::Pointer volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		typename TImageType::SpacingType spacing = volumeSkeleton->GetSpacing();

		BranchNode_t *portalNode = NULL;
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			if ( seedlist[0].toUInt() == m_branchNodeList.at(i)->RootPoint.x() &&
				seedlist[1].toUInt() == m_branchNodeList.at(i)->RootPoint.y() &&
				seedlist[2].toUInt() == m_branchNodeList.at(i)->RootPoint.z())
			{
				portalNode = m_branchNodeList.at(i); break;
			}
		}

		if ( portalNode == NULL)
		{
			qDebug() << "cannot locate portal point";
			return;
		}

		// remove extreme branches which near to portal node, two level
		// remove irregular branches of portal point
		for (int i = 0; i < portalNode->SubBranches.count(); i++)
		{
			if ( portalNode->SubBranchProperties[i]->extreme == true)
			{
				portalNode->SubBranchProperties[i]->useless = true;

				//m_displayModel->removeABranch(portalNode->id, i);
			}

			if ( portalNode->SubBranches[i]->count() < 50)
			{
				for ( int j = 0; j < portalNode->Children[i]->Children.count(); j++)
				{
					if ( portalNode->Children[i]->SubBranchProperties[j]->extreme == true || portalNode->Children[i]->maxDepth <=2)
					{
						portalNode->Children[i]->SubBranchProperties[j]->useless = true;

						//m_displayModel->removeABranch(portalNode->Children[i]->id, j);
					}
				}
			}
		}


		// remove extreme branches which less than 12 points (3 re-sampled point maybe)
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			for (int j = 0; j < m_branchNodeList[i]->SubBranches.count(); j++)
			{
				if ( m_branchNodeList[i]->SubBranches[j]->count()<=12 && m_branchNodeList[i]->SubBranchProperties[j]->extreme == true)
				{
					m_branchNodeList[i]->SubBranchProperties[j]->useless = true;

					//m_displayModel->removeABranch(m_branchNodeList[i]->id, j);
				}
			}

			m_branchNodeList[i]->group = 0;
		}



		int branchTypes = 0;
		int useful_portal_children_count = 0;
		QList<BranchNode_t*> nodes_group_left;
		QList<BranchNode_t*> nodes_group_right;

		// check type 3

		QList<BranchNode_t*> nodes_type3;

		for (int i = 0; i < portalNode->Children.count(); i++)
		{
			if ( isValidNode(portalNode->Children[i], 1, 3))
			{
				useful_portal_children_count++;
				nodes_type3.push_back(portalNode->Children[i]);
				qDebug() << " valid node : " << portalNode->Children[i]->RootPoint;
			}
		}

		if ( useful_portal_children_count == 3)
		{
			branchTypes = 3;
			qDebug() << "Chech type 3 [TRUE]" ;
			nodes_group_right.clear();
			nodes_group_left.clear();

			//grouping to left and right
			for (int i = 0; i < nodes_type3.count(); i++)
			{
				if ( nodes_type3[i]->RootPoint.x() < portalNode->RootPoint.x())
				{
					qDebug() << "right: " << nodes_type3[i]->RootPoint;
					nodes_group_right.push_back(nodes_type3[i]);
				}
				else
				{
					qDebug() << "left: " << nodes_type3[i]->RootPoint;
					nodes_group_left.push_back(nodes_type3[i]);
				}
			}

			if ( nodes_group_right.count() != 2)
			{
				qDebug() << "Type 3 Error: right branch has " << nodes_group_right.count() << " branches";
				return;
			}

			// process right
			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				nodes_group_right[i] = getNextValidNode(nodes_group_right[i], 2, 1);
				if ( nodes_group_right[i] == NULL)
				{
					qDebug() << "Type 3 Error, didn't get valid node";
					return;
				}
			}

			QList<BranchNode_t*> nodes_group_right_4case;
			for (int i = 0; i < nodes_group_right.count(); i++)
			{

				for (int j = 0; j < nodes_group_right[i]->Children.count(); j++)
				{
					//BranchNode_t* node = getNextValidNode(nodes_group_right[i]->Children[j], 2, 1);
					//if ( node)
					nodes_group_right_4case.push_back(nodes_group_right[i]->Children[j]);
				}
			}


			selectGroupingBranches<TImageType>(nodes_group_right_4case);

			return;
		}


		qDebug() << "Check type 3 : failed";


		// check type 4

		QList<BranchInfo_t*> nodes_level1;
		QList<BranchNode_t*> nodes_level2;
		nodes_group_left.clear();
		nodes_group_right.clear();

		for (int i = 0; i < portalNode->Children.count(); i++)
		{
			for (int j = 0; j < portalNode->Children[i]->Children.count(); j++)
			{
				nodes_level2.push_back(portalNode->Children[i]->Children[j]);
			}
		}

		for ( int i = 0; i < nodes_level2.count(); i++)
		{
			// 			if ( nodes_level2[i]->Children.count() < 1)
			// 				continue;

			if ( nodes_level2[i]->RootPoint.x() < portalNode->RootPoint.x())
			{
				qDebug() << "right" << nodes_level2[i]->RootPoint;
				nodes_group_right.push_back(nodes_level2[i]);
			}
			else
			{
				qDebug() << "left" << nodes_level2[i]->RootPoint;
				nodes_group_left.push_back(nodes_level2[i]);
			}
		}


		if (useful_portal_children_count == 2 && nodes_group_left.count() == 1 && nodes_group_right.count() == 3)
		{
			qDebug() << "It's possible branch type 4" << endl;
			BranchNode_t* goDeeper = NULL;

			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				if ( nodes_group_left[0]->Parent->RootPoint == nodes_group_right[i]->Parent->RootPoint)
				{
					branchTypes = 4;
					goDeeper = nodes_group_right[i];
					cout << "Branch type = " << branchTypes << endl;
				}
			}

			if ( branchTypes == 4 && goDeeper)
			{
				nodes_group_right.remove(goDeeper);

				qDebug() << "Go deeper children = " << goDeeper->RootPoint;

				//if ( goDeeper->Children.count() == 1)
				goDeeper = getNextValidNode(goDeeper, 2, 0);

				qDebug() << "Go deeper new = " << goDeeper->RootPoint;

				for (int i = 0; i < goDeeper->Children.count(); i++)
				{
					nodes_group_right.push_back(goDeeper->Children[i]);
				}

				if ( goDeeper->Children.count()>2)
				{
					BranchNode_t* maxDepthNode = goDeeper->Children[0];
					for (int i = 0; i < goDeeper->Children.count(); i++)
					{
						if ( goDeeper->Children[i]->maxDepth > maxDepthNode->maxDepth)
							maxDepthNode = goDeeper->Children[i];
					}

					maxDepthNode->group = 6;
					for (int i = 0; i < goDeeper->Children.count(); i++)
					{
						if ( goDeeper->Children[i] != maxDepthNode)
							goDeeper->Children[i]->group = 7;
					}
				}

				selectGroupingBranches<TImageType>(nodes_group_right);

				return;
			}
		}

		qDebug() << "Check type 4 : failed";


		// check type 2

		BranchNode_t* node_left;
		BranchNode_t* node_right;
		if ( useful_portal_children_count == 2)
		{
			//if ( portalNode->Children[0]->Children.count()> 0 && portalNode->Children[1]->Children.count()>0)
			{
				int children0_x_count = portalNode->Children[0]->RootPoint.x();
				int children1_x_count = portalNode->Children[1]->RootPoint.x();

				for (int i = 0; i < portalNode->Children[0]->Children.count(); i++)
					children0_x_count += portalNode->Children[0]->Children[i]->RootPoint.x();
				for (int i = 0; i < portalNode->Children[1]->Children.count(); i++)
					children1_x_count += portalNode->Children[1]->Children[i]->RootPoint.x();

				if ( children0_x_count < children1_x_count)
				{
					node_right = portalNode->Children[0];
					node_left = portalNode->Children[1];
				}
				else
				{
					node_right = portalNode->Children[1];
					node_left = portalNode->Children[0];
				}

				qDebug() << "left main node " << node_left->RootPoint << " x count" << children0_x_count;
				qDebug() << "right main node " << node_right->RootPoint << " x count" << children1_x_count;
			}
			// 			else
			// 			{
			// 				if ( portalNode->Children[0]->RootPoint.x() < portalNode->Children[1]->RootPoint.x())
			// 				{// right
			// 					node_right = portalNode->Children[0];
			// 					node_left = portalNode->Children[1];
			// 				}
			// 				else
			// 				{
			// 					node_right = portalNode->Children[1];
			// 					node_left = portalNode->Children[2];
			// 				}
			// 			}
		}

		BranchNode_t* node_right_new = NULL;
		node_right_new = getNextValidNode(node_right, 2, 2);

		if ( node_right_new == NULL)
		{
			qDebug() << "invalid node right, try less condition";

			node_right_new = getNextValidNode(node_right, 2, 1);

			if ( node_right_new == NULL)
			{
				qDebug() << "unknown type";
				return;
			}
		}

		qDebug() << "check type 2: valid next node " << node_right_new->RootPoint;
		qDebug() << "check type 2: has " << node_right_new->Children.count() << " children";
		nodes_group_left.clear();
		nodes_group_right.clear();

		if ( node_right_new->Children.count() >= 3)
		{
			qDebug() << "Possible type 2";
		}
		else if ( node_right_new->Children.count() == 2)
		{
			qDebug() << "Possible type 1";

			for (int i = 0; i < node_right_new->Children.count(); i++)
			{
				BranchNode_t* node = getNextValidNode(node_right_new->Children[i], 2, 1);

				for (int j = 0; j < node->Children.count();j++)
				{
					if ( node->SubBranchProperties[j]->useless == false)
					{
						nodes_group_right.push_back(node->Children[j]);
					}
				}
			}

			selectGroupingBranches<TImageType>(nodes_group_right);
		}
		else
		{
			qDebug() << "Unknown type";
			return;
		}












		return;

		// test validness of branches

		QList<BranchNode_t*>::iterator it_right = nodes_group_right.begin();
		//		QList<BranchNode_t*>::iterator it_left = nodes_group_left.begin();
		for (; it_right != nodes_group_right.end();)
		{
			if ( !isValidNode(*it_right, 1, 2))
			{
				qDebug() << "Invalid node [2 child,3 level] : " << (*it_right)->RootPoint;
				it_right = nodes_group_right.erase(it_right);
			}
			else
				it_right++;
		}

		// 		for (; it_left != nodes_group_left.end();)
		// 		{
		// 			if ( !isValidNode(*it_left, 2, 3))
		// 				it_left = nodes_group_right.erase(it_left)
		// 			else
		// 				it_right++;
		// 		}


		if ( nodes_group_right.count() == 1)
		{
			qDebug() << "NEED a valid node : " << nodes_group_right[0]->RootPoint;
			nodes_group_right[0] = getNextValidNode(nodes_group_right[0], 2, 3);

			if ( nodes_group_right[0] != NULL)
				qDebug() << "GOT a valid node : " << nodes_group_right[0]->RootPoint;
			else
				nodes_group_right.clear();
		}



		//relocate left branches
		if ( !branchTypes && nodes_group_left.count() == 1)
		{
			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				if ( nodes_group_right[i] == NULL)
					continue;

				if ( nodes_group_left[0]->Parent->RootPoint == nodes_group_right[i]->Parent->RootPoint)
				{
					branchTypes = 4;
					cout << "Branch type = " << branchTypes << endl;
				}
			}
		}

		if (!branchTypes)
		{
			if (nodes_group_right.count() == 3)
			{
				branchTypes = 2;
				cout << "Branch type = " << branchTypes << endl;
			}

			if (!branchTypes && nodes_group_right.count() == 4)
			{
				branchTypes = 1;
				cout << "Branch type = " << branchTypes << endl;
			}

			if (!branchTypes)
			{
				cout << "Branch type = unknown " << endl;
			}
		}


		int branch_color = 1;

		for ( int i = 0; i < nodes_group_right.count(); i++)
		{
			if ( nodes_group_right[i] == NULL)
				continue;

			QVector4D cvt_pnode = nodes_group_right[i]->Parent->RootPoint;
			cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
			cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
			cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

			QVector4D cvt_cnode =  nodes_group_right[i]->RootPoint;
			cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
			cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
			cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

			m_displayModel->setBranchClass(branch_color++,  cvt_pnode, cvt_cnode);
		}

		for ( int i = 0; i < nodes_group_left.count(); i++)
		{
			QVector4D cvt_pnode = nodes_group_left[i]->Parent->RootPoint;
			cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
			cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
			cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

			QVector4D cvt_cnode =  nodes_group_left[i]->RootPoint;
			cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
			cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
			cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

			m_displayModel->setBranchClass(branch_color++,  cvt_pnode, cvt_cnode);
		}


	}


	template<class TPixelType>
	void autoClassifyBranches()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		QStringList seedlist = autoLocatePortalPoint();
		if ( seedlist.size() != 4)
		{
			qDebug() << "[Error]: cannot find portal point";
			return;
		}

		typename TImageType::Pointer volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		typename TImageType::SpacingType spacing = volumeSkeleton->GetSpacing();

		BranchNode_t *portalNode = NULL;
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			if ( seedlist[0].toUInt() == m_branchNodeList.at(i)->RootPoint.x() &&
				 seedlist[1].toUInt() == m_branchNodeList.at(i)->RootPoint.y() &&
				 seedlist[2].toUInt() == m_branchNodeList.at(i)->RootPoint.z())
			{
				portalNode = m_branchNodeList.at(i); break;
			}
		}

		if ( portalNode == NULL)
		{
			qDebug() << "[Error]: cannot find portal point";
			return;
		}

		// remove extreme branches which near to portal node, two level
		// remove irregular branches of portal point
		for (int i = 0; i < portalNode->SubBranches.count(); i++)
		{
			if ( portalNode->SubBranchProperties[i]->extreme == true)
			{
				portalNode->SubBranchProperties[i]->useless = true;

				//m_displayModel->removeABranch(portalNode->id, i);
			}

			if ( portalNode->SubBranches[i]->count() < 50)
			{
				for ( int j = 0; j < portalNode->Children[i]->Children.count(); j++)
				{
					if ( portalNode->Children[i]->SubBranchProperties[j]->extreme == true || portalNode->Children[i]->maxDepth <=2)
					{
						portalNode->Children[i]->SubBranchProperties[j]->useless = true;

						//m_displayModel->removeABranch(portalNode->Children[i]->id, j);
					}
				}
			}
		}


		// remove extreme branches which less than 12 points (3 re-sampled point maybe)
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			for (int j = 0; j < m_branchNodeList[i]->SubBranches.count(); j++)
			{
				if ( m_branchNodeList[i]->SubBranches[j]->count()<=12 && m_branchNodeList[i]->SubBranchProperties[j]->extreme == true)
				{
					m_branchNodeList[i]->SubBranchProperties[j]->useless = true;

					m_displayModel->removeABranch(m_branchNodeList[i]->id, j);
				}
			}

			m_branchNodeList[i]->group = 0;
		}



		int branchTypes = 0;
		int useful_portal_children_count = 0;
		QList<BranchNode_t*> nodes_group_left;
		QList<BranchNode_t*> nodes_group_right;
		QList<BranchNode_t*> nodes_group_right_2parts;

		// check type 3

		QList<BranchNode_t*> nodes_type3;
		
		for (int i = 0; i < portalNode->Children.count(); i++)
		{
			if ( isValidNode(portalNode->Children[i], 1, 3))
			{
				useful_portal_children_count++;
				nodes_type3.push_back(portalNode->Children[i]);
				qDebug() << " valid node : " << portalNode->Children[i]->RootPoint;
			}
		}

		if ( useful_portal_children_count == 3)
		{
			branchTypes = 3;
			qDebug() << "Chech type 3 [TRUE]" ;
			nodes_group_right.clear();
			nodes_group_left.clear();

			//grouping to left and right
			for (int i = 0; i < nodes_type3.count(); i++)
			{
				if ( nodes_type3[i]->RootPoint.x() < portalNode->RootPoint.x())
				{
					qDebug() << "right: " << nodes_type3[i]->RootPoint;
					nodes_group_right.push_back(nodes_type3[i]);
				}
				else
				{
					qDebug() << "left: " << nodes_type3[i]->RootPoint;
					nodes_group_left.push_back(nodes_type3[i]);
				}
			}
			
			if ( nodes_group_right.count() != 2)
			{
				qDebug() << "Type 3 Error: right branch has " << nodes_group_right.count() << " branches";
				return;
			}

			// process right
			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				nodes_group_right[i] = getNextValidNode(nodes_group_right[i], 2, 1);
				if ( nodes_group_right[i] == NULL)
				{
					qDebug() << "Type 3 Error, didn't get valid node";
					return;
				}
			}

			/*
			QList<BranchNode_t*> nodes_group_right_4case;
			for (int i = 0; i < nodes_group_right.count(); i++)
			{

				for (int j = 0; j < nodes_group_right[i]->Children.count(); j++)
				{
					//BranchNode_t* node = getNextValidNode(nodes_group_right[i]->Children[j], 2, 1);
					//if ( node)
					nodes_group_right_4case.push_back(nodes_group_right[i]->Children[j]);
				}
			}*/

			
			selectGroupingBranches<TImageType>(nodes_group_right);

			return;
		}


		qDebug() << "Check type 3 : failed";


		// check type 4
		
		QList<BranchInfo_t*> nodes_level1;
		QList<BranchNode_t*> nodes_level2;
		nodes_group_left.clear();
		nodes_group_right.clear();

		for (int i = 0; i < portalNode->Children.count(); i++)
		{
			for (int j = 0; j < portalNode->Children[i]->Children.count(); j++)
			{
				nodes_level2.push_back(portalNode->Children[i]->Children[j]);
			}
		}

		for ( int i = 0; i < nodes_level2.count(); i++)
		{
// 			if ( nodes_level2[i]->Children.count() < 1)
// 				continue;

			if ( nodes_level2[i]->RootPoint.x() < portalNode->RootPoint.x())
			{
				qDebug() << "right" << nodes_level2[i]->RootPoint;
				nodes_group_right.push_back(nodes_level2[i]);
			}
			else
			{
				qDebug() << "left" << nodes_level2[i]->RootPoint;
				nodes_group_left.push_back(nodes_level2[i]);
			}
		}


		if (useful_portal_children_count == 2 && nodes_group_left.count() == 1 && nodes_group_right.count() == 3)
		{
			qDebug() << "It's possible branch type 4" << endl;
			BranchNode_t* goDeeper = NULL;

			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				if ( nodes_group_left[0]->Parent->RootPoint == nodes_group_right[i]->Parent->RootPoint)
				{
					branchTypes = 4;
					goDeeper = nodes_group_right[i];
					cout << "Branch type = " << branchTypes << endl;
					nodes_group_right_2parts.push_back(nodes_group_right[i]);
				}
				else
				{
					if ( nodes_group_right_2parts.find(nodes_group_right[i]->Parent) == nodes_group_right_2parts.end())
						nodes_group_right_2parts.push_back(nodes_group_right[i]->Parent);
				}
			}

			if ( nodes_group_right_2parts.count() == 2)
				selectGroupingBranches<TImageType>(nodes_group_right_2parts);

			/*
			if ( branchTypes == 4 && goDeeper)
			{
				nodes_group_right.remove(goDeeper);

				qDebug() << "Go deeper children = " << goDeeper->RootPoint;

				//if ( goDeeper->Children.count() == 1)
					goDeeper = getNextValidNode(goDeeper, 2, 0);

				qDebug() << "Go deeper new = " << goDeeper->RootPoint;

				for (int i = 0; i < goDeeper->Children.count(); i++)
				{
					nodes_group_right.push_back(goDeeper->Children[i]);
				}

				if ( goDeeper->Children.count()>2)
				{
					BranchNode_t* maxDepthNode = goDeeper->Children[0];
					for (int i = 0; i < goDeeper->Children.count(); i++)
					{
						if ( goDeeper->Children[i]->maxDepth > maxDepthNode->maxDepth)
							maxDepthNode = goDeeper->Children[i];
					}

					maxDepthNode->group = 6;
					for (int i = 0; i < goDeeper->Children.count(); i++)
					{
						if ( goDeeper->Children[i] != maxDepthNode)
							goDeeper->Children[i]->group = 7;
					}
				}

				selectGroupingBranches<TImageType>(nodes_group_right);

				return;
			}
			*/
		}

		qDebug() << "Check type 4 : failed";


		// check type 2 and 1

		BranchNode_t* node_left;
		BranchNode_t* node_right;
		if ( useful_portal_children_count == 2)
		{
			//if ( portalNode->Children[0]->Children.count()> 0 && portalNode->Children[1]->Children.count()>0)
			{
				int children0_x_count = portalNode->Children[0]->RootPoint.x();
				int children1_x_count = portalNode->Children[1]->RootPoint.x();

				for (int i = 0; i < portalNode->Children[0]->Children.count(); i++)
					children0_x_count += portalNode->Children[0]->Children[i]->RootPoint.x();
				for (int i = 0; i < portalNode->Children[1]->Children.count(); i++)
					children1_x_count += portalNode->Children[1]->Children[i]->RootPoint.x();

				if ( children0_x_count < children1_x_count)
				{
					node_right = portalNode->Children[0];
					node_left = portalNode->Children[1];
				}
				else
				{
					node_right = portalNode->Children[1];
					node_left = portalNode->Children[0];
				}

				qDebug() << "left main node " << node_left->RootPoint << " x count" << children0_x_count;
				qDebug() << "right main node " << node_right->RootPoint << " x count" << children1_x_count;
			}
// 			else
// 			{
// 				if ( portalNode->Children[0]->RootPoint.x() < portalNode->Children[1]->RootPoint.x())
// 				{// right
// 					node_right = portalNode->Children[0];
// 					node_left = portalNode->Children[1];
// 				}
// 				else
// 				{
// 					node_right = portalNode->Children[1];
// 					node_left = portalNode->Children[2];
// 				}
// 			}
		}

		BranchNode_t* node_right_new = NULL;
		node_right_new = getNextValidNode(node_right, 2, 2);

		if ( node_right_new == NULL)
		{
			qDebug() << "invalid node right, try less condition";

			node_right_new = getNextValidNode(node_right, 2, 1);

			if ( node_right_new == NULL)
			{
				qDebug() << "unknown type";
				return;
			}
		}

		qDebug() << "check type 2: valid next node " << node_right_new->RootPoint;
		qDebug() << "check type 2: has " << node_right_new->Children.count() << " children";
		nodes_group_left.clear();
		nodes_group_right.clear();

		if ( node_right_new->Children.count() >= 3)
		{
			qDebug() << "Possible type 2";

			BranchNode_t* node_ra = node_right_new->Children[0];

			for (int i = 0; i < node_right_new->Children.count(); i++)
			{
				if ( node_right_new->Children[i]->RootPoint.z() > node_ra->RootPoint.z())
					node_ra = node_right_new->Children[i];
			}


			node_ra->group = 8;
			
			nodes_group_right_2parts.clear();
			nodes_group_right_2parts.push_back(node_ra);
			for (int i = 0; i < node_right_new->Children.count(); i++)
			{
				if ( node_ra != node_right_new->Children[i])
				{
					node_right_new->Children[i]->group = 9;
					nodes_group_right_2parts.push_back(node_right_new->Children[i]);
				}
			}

			selectGroupingBranches<TImageType>(nodes_group_right_2parts);


		}
		else if ( node_right_new->Children.count() == 2)
		{
			qDebug() << "Possible type 1";

			nodes_group_right_2parts.clear();
			nodes_group_right_2parts.push_back(node_right_new->Children[0]);
			nodes_group_right_2parts.push_back(node_right_new->Children[1]);

			selectGroupingBranches<TImageType>(nodes_group_right_2parts);
		}
		else
		{
			qDebug() << "Unknown type";
			return;
		}












		return;

		// test validness of branches

		QList<BranchNode_t*>::iterator it_right = nodes_group_right.begin();
//		QList<BranchNode_t*>::iterator it_left = nodes_group_left.begin();
		for (; it_right != nodes_group_right.end();)
		{
			if ( !isValidNode(*it_right, 1, 2))
			{
				qDebug() << "Invalid node [2 child,3 level] : " << (*it_right)->RootPoint;
				it_right = nodes_group_right.erase(it_right);
			}
			else
				it_right++;
		}

// 		for (; it_left != nodes_group_left.end();)
// 		{
// 			if ( !isValidNode(*it_left, 2, 3))
// 				it_left = nodes_group_right.erase(it_left)
// 			else
// 				it_right++;
// 		}


		if ( nodes_group_right.count() == 1)
		{
			qDebug() << "NEED a valid node : " << nodes_group_right[0]->RootPoint;
			nodes_group_right[0] = getNextValidNode(nodes_group_right[0], 2, 3);
			
			if ( nodes_group_right[0] != NULL)
				qDebug() << "GOT a valid node : " << nodes_group_right[0]->RootPoint;
			else
				nodes_group_right.clear();
		}



		//relocate left branches
		if ( !branchTypes && nodes_group_left.count() == 1)
		{
			for (int i = 0; i < nodes_group_right.count(); i++)
			{
				if ( nodes_group_right[i] == NULL)
					continue;

				if ( nodes_group_left[0]->Parent->RootPoint == nodes_group_right[i]->Parent->RootPoint)
				{
					branchTypes = 4;
					cout << "Branch type = " << branchTypes << endl;
				}
			}
		}

		if (!branchTypes)
		{
			if (nodes_group_right.count() == 3)
			{
				branchTypes = 2;
				cout << "Branch type = " << branchTypes << endl;
			}
			
			if (!branchTypes && nodes_group_right.count() == 4)
			{
				branchTypes = 1;
				cout << "Branch type = " << branchTypes << endl;
			}
			
			if (!branchTypes)
			{
				cout << "Branch type = unknown " << endl;
			}
		}
		

		int branch_color = 1;

		for ( int i = 0; i < nodes_group_right.count(); i++)
		{
			if ( nodes_group_right[i] == NULL)
				continue;

			QVector4D cvt_pnode = nodes_group_right[i]->Parent->RootPoint;
			cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
			cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
			cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

			QVector4D cvt_cnode =  nodes_group_right[i]->RootPoint;
			cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
			cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
			cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

			m_displayModel->setBranchClass(branch_color++,  cvt_pnode, cvt_cnode);
		}

		for ( int i = 0; i < nodes_group_left.count(); i++)
		{
			QVector4D cvt_pnode = nodes_group_left[i]->Parent->RootPoint;
			cvt_pnode.setX(cvt_pnode.x()*m_currentSpacing[0]);
			cvt_pnode.setY(cvt_pnode.y()*m_currentSpacing[1]);
			cvt_pnode.setZ(cvt_pnode.z()*m_currentSpacing[2]);

			QVector4D cvt_cnode =  nodes_group_left[i]->RootPoint;
			cvt_cnode.setX(cvt_cnode.x()*m_currentSpacing[0]);
			cvt_cnode.setY(cvt_cnode.y()*m_currentSpacing[1]);
			cvt_cnode.setZ(cvt_cnode.z()*m_currentSpacing[2]);

			m_displayModel->setBranchClass(branch_color++,  cvt_pnode, cvt_cnode);
		}


	}


	template<class TPixelType>
	void autoDivideLeftLiver()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		typename TImageType::Pointer m_volumeCutting = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		typename TImageType::Pointer m_volumeLiver = ((VolumeSource<TImageType>*)m_sourceLiver)->getVolume();
		typename TImageType::SizeType size = m_volumeLiver->GetBufferedRegion().GetSize();

		vtkPoints *points = vtkPoints::New();
		typename TImageType::IndexType pixelIndex;
		for (unsigned int i = 0; i < size[2]; i++)
		{
			for (unsigned int j = 0; j< size[1]; j++ )
			{
				for (unsigned int k = 0; k < size[0]; k++)
				{
					pixelIndex[0] = k;
					pixelIndex[1] = j;
					pixelIndex[2] = i;

					if ( m_volumeCutting->GetPixel(pixelIndex) == 1)
					{
						points->InsertNextPoint((double)pixelIndex[0]*m_currentSpacing[0], 
							(double)pixelIndex[1]*m_currentSpacing[1], 
							(double)pixelIndex[2]*m_currentSpacing[2]);
					}
				}
			}
		}


		//----------------- testing

// 		SurfaceSource* liverSurface = volume2SmoothSurface(m_sourceLiver, QString("LiverSurface"), 1);
// 		liverSurface->sourceModel();
// 		m_sourceLiver->addChild(liverSurface);
// 
// 		vtkSelectPolyData *loop = vtkSelectPolyData::New();
// 		loop->SetInput(liverSurface->sdata());
// 		loop->SetLoop(points);
// 		loop->GenerateSelectionScalarsOn();
// 		loop->SetSelectionModeToSmallestRegion();
// 
// 		vtkClipPolyData *clip = vtkClipPolyData::New();
// 		clip->SetInputConnection(loop->GetOutputPort());
// 
// 		vtkPolyDataMapper *clipMapper = vtkPolyDataMapper::New();
// 		clipMapper->SetInputConnection(clip->GetOutputPort());
// 
// 		vtkLODActor *clipActor = vtkLODActor::New();
// 		clipActor->SetMapper(clipMapper);
// 
// 		m_displayModel->addActor(clipActor);
// 
// 		//------------------ end testing
// 
// 		return;








		QStringList seedlist = autoLocatePortalPoint();
		if ( seedlist.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return;
		}
		else
		{
			qDebug() << "choose seed point: " << seedlist;
		}

		typename TImageType::Pointer volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		//typename TImageType::Pointer volumeLiver = ((VolumeSource<TImageType>*)m_sourceLiver)->getVolume();
		typename TImageType::SpacingType spacing = volumeSkeleton->GetSpacing();


		BranchNode_t *portalNode = NULL;

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			if ( seedlist[0].toUInt() == m_branchNodeList.at(i)->RootPoint.x() &&
				seedlist[1].toUInt() == m_branchNodeList.at(i)->RootPoint.y() &&
				seedlist[2].toUInt() == m_branchNodeList.at(i)->RootPoint.z())
			{
				portalNode = m_branchNodeList.at(i);
				break;
			}
		}

		if ( portalNode == NULL)
		{
			qDebug() << "cannot locate portal point";
			return;
		}

		QQueue<BranchNode_t*> nodeQueue;
		nodeQueue.push_back(portalNode);
		BranchNode_t* portalBranchNode = NULL;
		while(1)
		{
			bool bGotPortalPoint = true;
			if ( nodeQueue.count() < 1)
				break;

			BranchNode_t* currnode = nodeQueue.front();
			nodeQueue.pop_front();
			QList<int> child_dist_list;

			for (int i = 0; i < currnode->SubBranches.count(); i++)
			{
				QList<QVector4D> *branch = currnode->SubBranches[i];
				BranchProperty_t *branchProperty = currnode->SubBranchProperties[i];

				if ( branchProperty->useless || branchProperty->extreme)
					continue;

				QVector4D v4 = branch->at(branch->count()-1);
				int dist = sqrt((currnode->RootPoint.x()-v4.x())*(currnode->RootPoint.x()-v4.x())
					+ (currnode->RootPoint.y()-v4.y())*(currnode->RootPoint.y()-v4.y())
					+ (currnode->RootPoint.z()-v4.z())*(currnode->RootPoint.z()-v4.z()));

				qDebug() << "distance to child " << i << " " << dist << endl;
				child_dist_list.push_back(dist);
			}

			if ( child_dist_list.count() < 2)
			{
				bGotPortalPoint = false;
			}
			else
			{
				for (int i = 0; i < child_dist_list.count(); i++)
				{
					if ( child_dist_list[i] < 5)
					{
						bGotPortalPoint = false;
						break;
					}
				}
			}

			if ( bGotPortalPoint == true)
			{
				QString str = QString("%1 %2 %3 %4").arg(currnode->RootPoint.x()).arg(currnode->RootPoint.y()).arg(currnode->RootPoint.z()).arg(currnode->RootPoint.w());
				m_portalPoint = str.split(" ");
				portalBranchNode = currnode;
				break;
			}
			else 
			{
				for (int j = 0; j < currnode->Children.count();j++)
					nodeQueue.push_back(currnode->Children.at(j));
			}
		}

		m_displayModel->addSphere(m_portalPoint[0].toDouble()*spacing[0], 
			m_portalPoint[1].toDouble()*spacing[1], 
			m_portalPoint[2].toDouble()*spacing[2], m_portalPoint[3].toDouble()*0.5);	
		qDebug() << "Got portal point " << m_portalPoint;


		//find out left point

		BranchNode_t* leftNode1 = portalBranchNode->Children[0];
		BranchNode_t* leftNode1_1 = NULL;

		for (int i = 0; i < portalBranchNode->Children.count();i++)
		{
			if(leftNode1->RootPoint.x() < portalBranchNode->Children[i]->RootPoint.x())
				leftNode1 = portalBranchNode->Children[i];
		}

		m_displayModel->addSphere(leftNode1->RootPoint.x()*spacing[0], 
			leftNode1->RootPoint.y()*spacing[1], 
			leftNode1->RootPoint.z()*spacing[2], leftNode1->RootPoint.w()*0.5);

		if ( leftNode1->Children.count() < 2)
		{
			qDebug() << "left node is not valid";
			return;
		}

		qDebug() << "left node has " << leftNode1->Children.count() << " children";

		
		QList<BranchNode_t*> validChildren;

		for (int i = 0; i < leftNode1->Children.count(); i++)
		{
			bool bExist = false;
			for(int j = 0; j < validChildren.count(); j++)
			{
				if ( leftNode1->Children.at(i)->RootPoint == validChildren.at(j)->RootPoint)
				{
					bExist = true;
					break;
				}
			}

			if ( !bExist)
			{
				validChildren.push_back(leftNode1->Children.at(i));
				qDebug() << leftNode1->Children[i]->RootPoint;
			}
		}


		qDebug() << "it has " << validChildren.count() << " valid children ";
		if ( validChildren.count()!=2)
		{
			qDebug() << " exit because of the children number" ;
			return;
		}

		for (int i = 0; i < validChildren.count(); i++)
			qDebug() << "  " << validChildren.at(i)->RootPoint;

		qDebug() << "left node has " << leftNode1->SubBranches.count() << " branches";
		for (int i = 0; i < leftNode1->Children.count(); i++)
		{
			qDebug() << "  useless: " << leftNode1->SubBranchProperties.at(i)->useless ;
			qDebug() << "  extreme: " << leftNode1->SubBranchProperties.at(i)->extreme ;
			QList<QVector4D>* subbranch = leftNode1->SubBranches.at(i);
			qDebug() << "  point  : " << subbranch->at(subbranch->count()-1);
		}

		QVector3D direction = (leftNode1->RootPoint - portalNode->RootPoint).toVector3D();
		QVector3D p1 = validChildren[0]->RootPoint.toVector3D();
		QVector3D p2 = validChildren[1]->RootPoint.toVector3D();

		double dist1 = p1.distanceToLine(leftNode1->RootPoint.toVector3D(), direction.normalized());
		double dist2 = p2.distanceToLine(leftNode1->RootPoint.toVector3D(), direction.normalized());

		qDebug() << "distance 1 = " << dist1 << endl;
		qDebug() << "distance 2 = " << dist2 << endl;

		QList<QVector4D> *branch_div1 = NULL;
		BranchProperty_t *branchProperty_div1 = NULL;

		//if ( dist1 > dist2)
		if ( validChildren[0]->RootPoint.y() < validChildren[1]->RootPoint.y())
		{
			leftNode1_1 = validChildren[0];	//then branch[1] is the DIV_1

			for (int i = 0; i < leftNode1->SubBranches.count(); i++)
			{
				QList<QVector4D>* branch = leftNode1->SubBranches.at(i);
				if ( validChildren[1]->RootPoint.x() == (branch->at(branch->count()-1)).x() &&
					 validChildren[1]->RootPoint.y() == (branch->at(branch->count()-1)).y() &&
					 validChildren[1]->RootPoint.z() == (branch->at(branch->count()-1)).z())
				{
					branch_div1 = leftNode1->SubBranches.at(i);
					branchProperty_div1 = leftNode1->SubBranchProperties.at(i);
					qDebug() << "got branch";
				}
			}
			
			
			m_displayModel->addSphere(validChildren[0]->RootPoint.x()*spacing[0], 
									  validChildren[0]->RootPoint.y()*spacing[1], 
									  validChildren[0]->RootPoint.z()*spacing[2], 
									  validChildren[0]->RootPoint.w()*0.5);
		}
		else
		{
			leftNode1_1 = validChildren[1];
			for (int i = 0; i < leftNode1->SubBranches.count(); i++)
			{
				QList<QVector4D>* branch = leftNode1->SubBranches.at(i);
				if ( validChildren[0]->RootPoint.x() == (branch->at(branch->count()-1)).x() &&
					validChildren[0]->RootPoint.y() == (branch->at(branch->count()-1)).y() &&
					validChildren[0]->RootPoint.z() == (branch->at(branch->count()-1)).z())
				{
					branch_div1 = leftNode1->SubBranches.at(i);
					branchProperty_div1 = leftNode1->SubBranchProperties.at(i);
					qDebug() << "got branch";
				}
			}


			m_displayModel->addSphere(validChildren[1]->RootPoint.x()*spacing[0], 
									  validChildren[1]->RootPoint.y()*spacing[1], 
									  validChildren[1]->RootPoint.z()*spacing[2], 
									  validChildren[1]->RootPoint.w()*0.5);
		}

		if (m_ui.m_btnWithLiver->isChecked() == false)
			return;

		//for division 1

		if ( !branch_div1)
			return;

		QList<QVector4D> ql = convertVolumeBranch2DisplayBranch<TImageType>(branch_div1, branchProperty_div1, spacing);
		m_displayModel->setSelectABranch(ql);

		divideBySelectedBranch<TImageType>(0);
		m_displayModel->clearSelectedBranches();

		// for division 2 and 3
		QList<QList<QVector4D>*> branches_div2; QList<BranchProperty_t*> branches_properties_div2;
		QList<QList<QVector4D>*> branches_div3; QList<BranchProperty_t*> branches_properties_div3;
		
		for (int i = 0; i < leftNode1_1->Children.count(); i++)
		{
			for (int j = 0; j < leftNode1_1->SubBranches.count(); j++)
			{
				QList<QVector4D>* branch = leftNode1_1->SubBranches.at(j);
				BranchProperty_t* branch_properties = leftNode1_1->SubBranchProperties.at(j);

				if ( leftNode1_1->Children[i]->RootPoint.x() == (branch->at(branch->count()-1)).x() &&
					leftNode1_1->Children[i]->RootPoint.y() == (branch->at(branch->count()-1)).y() &&
					leftNode1_1->Children[i]->RootPoint.z() == (branch->at(branch->count()-1)).z() )
				{
					if ( leftNode1_1->RootPoint.x() > leftNode1_1->Children[i]->RootPoint.x())
					{
						branches_div2.push_back(branch);	
						branches_properties_div2.push_back(branch_properties);
					}
					else
					{
						branches_div3.push_back(branch);
						branches_properties_div3.push_back(branch_properties);	
					}
				}
			}
		}



		qDebug() << "division 2 include " << branches_div2.count() << " branches";
		qDebug() << "division 3 include " << branches_div3.count() << " branches";

		for (int i = 0; i < branches_div2.count(); i++)
		{
			QList<QVector4D> ql = convertVolumeBranch2DisplayBranch<TImageType>(branches_div2.at(i), branches_properties_div2.at(i), spacing);
			m_displayModel->setSelectABranch(ql, 1);
		}

		divideBySelectedBranch<TImageType>(1);
		m_displayModel->clearSelectedBranches();

		for (int i = 0; i < branches_div3.count(); i++)
		{
			QList<QVector4D> ql = convertVolumeBranch2DisplayBranch<TImageType>(branches_div3.at(i), branches_properties_div3.at(i), spacing);
			m_displayModel->setSelectABranch(ql, 2);
		}

		divideBySelectedBranch<TImageType>(2);
		m_displayModel->clearSelectedBranches();

		//auto divide liver


	}


	template<class TImageType>
	void divideBySelectedBranch(int color)
	{
		typedef itk::Image<unsigned char, 3> TImageUChar3D;
		typedef itk::CastImageFilter<TImageType, TImageUChar3D> CastImageFilterShort2UChar;
		typename TImageType::Pointer orig_volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		typename TImageType::Pointer orig_volumeLiver = ((VolumeSource<TImageType>*)m_sourceLiver)->getVolume();
		typename TImageType::SpacingType spacing = orig_volumeLiver->GetSpacing();

		typename TImageUChar3D::Pointer volumeSkeleton;
		typename TImageUChar3D::Pointer volumeLiver;

		CastImageFilterShort2UChar::Pointer castSk = CastImageFilterShort2UChar::New();
		castSk->SetInput(orig_volumeSkeleton);
		castSk->Update();
		volumeSkeleton = castSk->GetOutput();

		CastImageFilterShort2UChar::Pointer castLi = CastImageFilterShort2UChar::New();
		castLi->SetInput(orig_volumeLiver);
		castLi->Update();
		volumeLiver = castLi->GetOutput();

		QList<QVector4D> seeds;
		QList<QVector4D> directions;
		QList<QList<QVector4D>> selectedBranches = m_displayModel->getSelectedBranch();	//get selected color

		for (int i = 0; i < selectedBranches.count(); i++)
		{
			QVector4D qseed = selectedBranches.at(i).at(2);
			qseed.setX(qseed.x()/spacing[0]);
			qseed.setY(qseed.y()/spacing[1]);
			qseed.setZ(qseed.z()/spacing[2]);
			seeds.push_back(qseed);

			QVector4D qseedDirection = selectedBranches.at(i).at(3);
			qseedDirection.setX(qseedDirection.x()/spacing[0]);
			qseedDirection.setY(qseedDirection.y()/spacing[1]);
			qseedDirection.setZ(qseedDirection.z()/spacing[2]);
			directions.push_back(qseedDirection);
		}

		divideLiverBySeedWithDirections<TImageUChar3D>(1,seeds, directions, volumeSkeleton, volumeLiver, color);
	}




	template<class TImageType>
	QList<QVector4D> convertVolumeBranch2DisplayBranch(QList<QVector4D>* branch, BranchProperty_t* branchProperty, typename TImageType::SpacingType spacing)
	{
		QList<QVector4D> ql;
		for (int i = 0; i < branch->count(); i++)
		{
			QVector4D vec = QVector4D((*branch)[i].x()*spacing[0],	(*branch)[i].y()*spacing[1], (*branch)[i].z()*spacing[2], (*branch)[i].w());

 			if ( i == branch->count()-1 && branchProperty->extreme == true)
 				vec.setW(0);

			ql.append(vec);
		}

		return ql;
	}


	template<class TPixelType>
	void autoDivideRightLiver()
	{

	}



	template<class TPixelType>
	bool rebuildVesselStructure()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		if ( m_sourceSkeleton == NULL)
		{
			qDebug() << "[warning] please choose skeleton volume";
			return false;
		}

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();//((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename TImageType::SizeType inputImageSize = m_inputImage->GetBufferedRegion().GetSize();
		typename TImageType::SpacingType spacing = m_inputImage->GetSpacing();
		cout << "spacing " << spacing[0] << ", " << spacing[1] << ", " << spacing[2] << endl;

		if ( m_ui.m_cbSeedList->count() == 0)
		{
			cout << "please load node points" << endl;
			return false;
		}

		QStringList seedlist = m_ui.m_cbSeedList->currentText().split(" ");//autoLocatePortalPoint<TPixelType>();
		if ( seedlist.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return false;
		}
		else
		{
			qDebug() << "choose seed point: " << seedlist;
		}

		int startingSeed[3];
		double startingRadius = seedlist[3].toDouble();

		startingSeed[0] = seedlist[0].toUInt();
		startingSeed[1] = seedlist[1].toUInt();
		startingSeed[2] = seedlist[2].toUInt();


		typename ImageDuplicator::Pointer dup0 = ImageDuplicator::New();
		dup0->SetInputImage(m_inputImage);
		dup0->Update();
		dup0->GetOutput()->FillBuffer(0);
		typename TImageType::Pointer m_outputImage = dup0->GetOutput();
		

		QVector4D startPointV4;
		startPointV4.setX(int(seedlist[0].toUInt()));
		startPointV4.setY(int(seedlist[1].toUInt()));
		startPointV4.setZ(int(seedlist[2].toUInt()));
		startPointV4.setW(startingRadius);

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			delete m_branchNodeList.at(i);
		}
		m_branchNodeList.clear();


		QQueue<QVector4D>	rootStack;		// a group stack is to find out a root and its children
		rootStack.enqueue(startPointV4);

		while(1)
		{
			if ( rootStack.count()==0)
				break;

			QVector4D rootPoint = rootStack.dequeue();		// the first point in a group is the root point
			BranchNode_t *branchPoint_t = new BranchNode_t;	// store the group
			branchPoint_t->removed = false;
			branchPoint_t->RootPoint = rootPoint;
			m_branchNodeList.push_back(branchPoint_t);
			QList<QVector4D> branches;

			// find out how many branches
			int radius = 2;
			for ( unsigned int l = rootPoint.x()-radius; l <= rootPoint.x()+radius; l++)
			{
				for ( unsigned int m = rootPoint.y()-radius; m <= rootPoint.y()+radius; m++)
				{
					for ( unsigned int n = rootPoint.z()-radius; n <= rootPoint.z()+radius; n++)
					{
						typename TImageType::IndexType neighbourPointIndex;
						neighbourPointIndex[0] = l;
						neighbourPointIndex[1] = m;
						neighbourPointIndex[2] = n;

						if ( l == rootPoint.x() && m == rootPoint.y() && n==rootPoint.z())
						{
							continue;
						}

						//only use the surface to estimate branches
						if ( std::labs(l-rootPoint.x()) == radius || 
							 std::labs(m-rootPoint.y()) == radius || 
							 std::labs(n-rootPoint.z())== radius)
						{
							if ( m_inputImage->GetPixel(neighbourPointIndex) == 1)	// has child
							{
								QVector4D p; p.setX(l); p.setY(m); p.setZ(n);

								branches.push_back(p);
								//rootStack.enqueue(p);
							}
						}
					}
				}
			}

			//cout << "branch number " << branches.count() << endl;

			// fill the output image
			for ( unsigned int l = rootPoint.x()-radius; l <= rootPoint.x()+radius; l++)
			{
				for ( unsigned int m = rootPoint.y()-radius; m <= rootPoint.y()+radius; m++ )
				{
					for ( unsigned int n = rootPoint.z()-radius; n <= rootPoint.z()+radius; n++)
					{
						typename TImageType::IndexType neighbourPointIndex;
						neighbourPointIndex[0] = l;
						neighbourPointIndex[1] = m;
						neighbourPointIndex[2] = n;

						if ( m_inputImage->GetPixel(neighbourPointIndex) == 1)
							m_outputImage->SetPixel(neighbourPointIndex, 1);
					}
				}
			}
			


			// for each child, growing until a 3-branch i < 1; i++ )//
			for (int i = 0; i < branches.count(); i++)
			{
				QStack<QVector4D> tempStack;
				tempStack.push(branches.at(i));
				QList<QVector4D>* branch = new QList<QVector4D>();
				BranchProperty_t* branchproperty = new BranchProperty_t;
				branchPoint_t->SubBranches.push_back(branch);
				branchPoint_t->SubBranchProperties.push_back(branchproperty);
				branch->push_back(rootPoint);

				while(1)
				{
					if ( tempStack.count() == 0)
						break;

					QVector4D centerPoint = tempStack.pop();
					typename TImageType::IndexType centerIndex;
					centerIndex[0] = centerPoint.x();
					centerIndex[1] = centerPoint.y();
					centerIndex[2] = centerPoint.z();

					//cout << "centerIndex: " << centerIndex[0] << ", " << centerIndex[1] << ", "
					//	<< centerIndex[2] << endl;

					int adjacentNum = 0;
					for ( unsigned int l = centerIndex[0]-1; l <= centerIndex[0]+1; l++)
					{
						for ( unsigned int m =  centerIndex[1]-1; m <=  centerIndex[1]+1; m++ )
						{
							for ( unsigned int n =  centerIndex[2]-1; n <=  centerIndex[2]+1; n++)
							{
								typename TImageType::IndexType neighbourIndex;
								neighbourIndex[0] = l;
								neighbourIndex[1] = m;
								neighbourIndex[2] = n;

								if ( l == centerIndex[0] && m == centerIndex[1] && n == centerIndex[2])
									continue;

								if ( m_inputImage->GetPixel(neighbourIndex) == 1)	// has child
								{
									if (m_outputImage->GetPixel(neighbourIndex) == 0)
									{
										QVector4D p(l,m,n,0);
										tempStack.push(p);
									}
									adjacentNum++;
								}
							}
						}
					}


					// add the radius information of a node
					for (int i = 0; i < m_ui.m_cbSeedList->count(); i++)
					{
						QStringList sl = (m_ui.m_cbSeedList->itemText(i)).split(" ");
						if (std::labs(sl[0].toUInt()-centerPoint.x()) <= 2 && 
							std::labs(sl[1].toUInt()-centerPoint.y()) <= 2 && 
							std::labs(sl[2].toUInt()-centerPoint.z()) <= 2)
						{
							centerPoint.setW(sl[3].toDouble());
							break;
						}
					}

					branch->push_back(centerPoint);
					//branchPoint_t->PointsList.push_back(currentPoint);
					m_outputImage->SetPixel(centerIndex, 1);

					if ( adjacentNum >= 3||adjacentNum == 1) // reach to another branch point
					{
						if ( adjacentNum == 1)
							branchproperty->extreme = true;
						else
							branchproperty->extreme = false;

						//ignore the root which appeared
						bool ignore = false;
						for (int j = 0; j < m_branchNodeList.count(); j++)
						{
							if ( m_branchNodeList.at(j)->RootPoint == centerPoint)
							{
								ignore = true;
								break;
							}
						}

						if ( !ignore)
						{
							rootStack.enqueue(centerPoint);


							break; //break while
						}//if
					}//if
				} // while
			}//for	
		}//while




		//check branches to delete
		int combiningThreshold = 4;

		for ( int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* branchnode = m_branchNodeList[i];

			for (int j = 0; j < branchnode->SubBranches.count(); j++)
			{
				QList<QVector4D>* subbranchpoints = branchnode->SubBranches[j];
				BranchProperty_t* branchproperty = branchnode->SubBranchProperties[j];
				branchproperty->useless = false;

				if ( subbranchpoints->count() < combiningThreshold)
				{
					branchproperty->useless = true;
				}
			}
		}


		//combine branches
		for ( int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* currentNode = m_branchNodeList[i];
			currentNode->removed = true;
			QList<QList<QVector4D>*> newChildBranches;
			QList<BranchProperty_t*> newChildProperty;

			for (int j = 0; j < currentNode->SubBranches.count(); j++)
			{
				QList<QVector4D>* branch = currentNode->SubBranches[j];
				QVector4D subbranch_tearnode = (*branch)[branch->count()-1];	// get the last node of this branch

				BranchProperty_t* branchproperty = currentNode->SubBranchProperties[j];

				if ( branchproperty->useless == false)
					continue;

				//search the node who are using the tearnode as ROOT, then remove the node
				for (int k = 0; k < m_branchNodeList.count(); k++)
				{
					BranchNode_t* bn = m_branchNodeList[k];
					QVector4D rootPoint = bn->RootPoint;

					if ( subbranch_tearnode == currentNode->RootPoint)		// don't move itself
					{
						//cout << "it is myself" << endl;
						continue;
					}

					if ( std::labs(subbranch_tearnode.x() - rootPoint.x()) <= 1 && 
						std::labs(subbranch_tearnode.y() - rootPoint.y()) <= 1 && 
						std::labs(subbranch_tearnode.z() - rootPoint.z()) <= 1)
					{
						cout << "find out a root point which is the child root of this branch" << endl;

						//change all root of new branches to this root
						for (int m = 0; m < bn->SubBranches.count(); m++)
						{
							QList<QVector4D>* bnsubbranch = bn->SubBranches[m];
							(*bnsubbranch)[0] = currentNode->RootPoint;
							newChildBranches.push_back(bnsubbranch);
							newChildProperty.push_back(bn->SubBranchProperties[m]);
						}

						bn->RootPoint = currentNode->RootPoint;
						//bn->removed = true;
					}
				}
			}

			for (int i = 0; i < newChildBranches.count(); i++)
			{
				currentNode->SubBranches.push_back(newChildBranches[i]);
				currentNode->SubBranchProperties.push_back(newChildProperty[i]);
			}
		}


		//construct hierarchical

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* currentNode = m_branchNodeList[i];
			currentNode->Parent = NULL;

			//qDebug() << "parent: " << currentNode->RootPoint;
			//qDebug() << "braNum: " << currentNode->SubBranches.count();

			for (int j = 0; j < currentNode->SubBranches.count(); j++)
			{
				QList<QVector4D>* branch = currentNode->SubBranches[j];
				QVector4D subbranch_tearnode = (*branch)[branch->count()-1];	// get the last node of this branch

				BranchProperty_t* branchproperty = currentNode->SubBranchProperties[j];

				if ( branchproperty->useless == true)
					continue;

				for (int k = 0; k < m_branchNodeList.count(); k++)
				{
					BranchNode_t *hitNode = m_branchNodeList[k];
					//if ( hitNode->removed)
					//	continue;

					if ( hitNode != currentNode)
					{
						if ( std::labs(subbranch_tearnode.x() - hitNode->RootPoint.x()) <= 1 && 
							std::labs(subbranch_tearnode.y() - hitNode->RootPoint.y()) <= 1 && 
							std::labs(subbranch_tearnode.z() - hitNode->RootPoint.z()) <= 1)
						{
							bool bExist = false;
							for (int m = 0; m < currentNode->Children.count(); m++)
							{
								if (currentNode->Children.at(m)->RootPoint == hitNode->RootPoint)
									bExist = true;
							}

							if ( bExist == false)
							{
								//qDebug() << "child : " << hitNode->RootPoint ;

								currentNode->Children.push_back(hitNode);
							}
						}
					}
				}
			}
		}



		//remove duplicated nodes

		QList<BranchNode_t*> newBranchList;
		for(int i = 0; i < m_branchNodeList.count();i++)
		{
			newBranchList.push_back(m_branchNodeList.at(i));
		}

		m_branchNodeList.clear();
		for(int i = 0; i < newBranchList.count(); i++)
		{
			bool bExist = false;
			for(int j=0; j < m_branchNodeList.count(); j++)
			{
				if ( m_branchNodeList.at(j)->RootPoint == newBranchList.at(i)->RootPoint)
				{
					bExist = true;	// but need add children info of this node to 

					//complete child nodes
					for (int m = 0; m < newBranchList.at(i)->Children.count(); m++)
					{
						bool bChildExist = false;
						BranchNode_t* node1 = newBranchList.at(i)->Children.at(m);

						for (int n = 0; n < m_branchNodeList.at(j)->Children.count(); n++)
						{
							BranchNode_t* node2 = m_branchNodeList.at(j)->Children.at(n);

							if ( node2->RootPoint == node1->RootPoint)
							{
								bChildExist = true;
								break;
							}
						}

						if ( !bChildExist)
							m_branchNodeList.at(j)->Children.push_back(newBranchList.at(i)->Children.at(m));
					}

					//complete child branches

					for (int m = 0; m < newBranchList.at(i)->SubBranches.count(); m++)
					{
						bool bBranchExist = false;
						QList<QVector4D>* branch1 = newBranchList.at(i)->SubBranches.at(m);

						for (int n = 0; n < m_branchNodeList.at(j)->SubBranches.count(); n++)
						{
							QList<QVector4D>* branch2 = m_branchNodeList.at(j)->SubBranches.at(n);

							if ( branch2->count() == branch1->count() && 
								branch2->at(branch2->count()-1) == branch1->at(branch1->count()-1))
							{

								bBranchExist = true;
								break;
							}	  
						}

						if ( !bBranchExist)
						{
							m_branchNodeList.at(j)->SubBranches.push_back(newBranchList.at(i)->SubBranches.at(m));
							m_branchNodeList.at(j)->SubBranchProperties.push_back(newBranchList.at(i)->SubBranchProperties.at(m));
						}
					}

					//printf("m_node_child_num %d, new_node_child_num %d\n", m_branchNodeList.at(j)->Children.count(), newBranchList.at(i)->Children.count());
					break;
				}
			}

			if ( !bExist)
				m_branchNodeList.push_back(newBranchList.at(i));
		}


		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* currentNode = m_branchNodeList[i];

			currentNode->id = i;

			for (int j = 0; j < currentNode->Children.count(); j++)
			{
				currentNode->Children.at(j)->Parent = currentNode;
			}
		}


		m_ui.m_cbSeedList->clear();

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* currentNode = m_branchNodeList[i];

			QString str = QString("%1 %2 %3 %4").
				arg(currentNode->RootPoint.x()).
				arg(currentNode->RootPoint.y()).
				arg(currentNode->RootPoint.z()).
				arg(currentNode->RootPoint.w());
				
			m_ui.m_cbSeedList->addItem(str);

			cout << currentNode->id << "\t";
			if ( currentNode->Parent)
				cout << currentNode->Parent->id << "\t";
			else
				cout << "-1\t" ;
			cout << currentNode->RootPoint.x() << "\t" 
				 << currentNode->RootPoint.y() << "\t" 
				 << currentNode->RootPoint.z() << "\t"
				 << currentNode->RootPoint.w() << "\t";

			BranchNode_t* parentNode = NULL;
			if ( currentNode->Parent)
			{
				for(int j = 0; j < m_branchNodeList.count(); j++)
				{
					if (currentNode->Parent->id == m_branchNodeList.at(j)->id)
					{
						parentNode = m_branchNodeList.at(j);
						break;
					}
				}
			}

			int branchlength = 0;
			if ( parentNode)
			{
				for (int m = 0; m < parentNode->SubBranches.count(); m++)
				{
					QList<QVector4D>* subbranch = parentNode->SubBranches.at(m);
					if (subbranch->at(subbranch->count()-1) == currentNode->RootPoint)
					{
						branchlength = subbranch->count();
						break;
					}
				}
			}

			cout << "\t" << branchlength << "\t" << endl;
		}


		// remove useless branches --------------------

		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* bp = m_branchNodeList[i];
			QList<QList<QVector4D>*> tmpSubbranches(bp->SubBranches);
			QList<BranchProperty_t*> tmpSubbranchProperties(bp->SubBranchProperties);
			bp->SubBranches.clear();
			bp->SubBranchProperties.clear();

			for (int c = 0; c < bp->Children.size(); c++)
			{
				for (int j = 0; j < tmpSubbranches.count(); j++)
				{
					QList<QVector4D> *branch = tmpSubbranches[j];
					BranchProperty_t *branchProperty = tmpSubbranchProperties[j];

					if ( branchProperty->useless == true)
						continue;

					qDebug() << branch->at(branch->size()-1) << ", " << bp->Children[c]->RootPoint ;

					if (branch->at(branch->size()-1) == bp->Children[c]->RootPoint)
					{
						bp->SubBranches.push_back(branch);
						bp->SubBranchProperties.push_back(branchProperty);
					}
				}
			}
		}

		// ------------------------------------------

		// calculate depth
		
		//display
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			BranchNode_t* bp = m_branchNodeList[i];

			for (int j = 0; j < bp->SubBranches.count(); j++)
			{
				QList<QVector4D> *branch = bp->SubBranches[j];
				BranchProperty_t *branchProperty = bp->SubBranchProperties[j];

				if ( branchProperty->useless)
					continue;

				QList<QVector4D> ql;
				for (int i = 0; i < branch->count(); i++)
				{
					QVector4D vec = QVector4D((*branch)[i].x()*spacing[0],
						(*branch)[i].y()*spacing[1],
						(*branch)[i].z()*spacing[2], 
						(*branch)[i].w());

					if ( i == branch->count()-1 && branchProperty->extreme == true)
						vec.setW(0);

					ql.append(vec);
				}
				m_displayModel->addBranch(ql, bp->id, j );
			}
		}

		m_displayModel->resetCamera();
		m_displayModel->updateActive3DScene();

		calculateMaxDepth();
		calculateDepthFromPortal();
		slotLocateHepaticVein();

		if (m_ui.m_cbRemoveHV->isChecked())
			slotRemoveHepatic();

		return true;
	}


	void calculateDepthFromPortal()
	{
		QStringList seedlist = autoLocatePortalPoint();
		if ( seedlist.size() != 4)
		{
			qDebug() << "[Error]: cannot find portal point";
			return;
		}

		BranchNode_t *portalNode = NULL;
		for (int i = 0; i < m_branchNodeList.count(); i++)
		{
			if ( seedlist[0].toUInt() == m_branchNodeList.at(i)->RootPoint.x() &&
				seedlist[1].toUInt() == m_branchNodeList.at(i)->RootPoint.y() &&
				seedlist[2].toUInt() == m_branchNodeList.at(i)->RootPoint.z())
			{
				portalNode = m_branchNodeList.at(i); break;
			}
		}

		if ( portalNode == NULL)
		{
			qDebug() << "[Error]: cannot find portal point";
			return;
		}

		int currentLevel = 0;
		QList<BranchNode_t*> node_list;
		portalNode->depthFromPortal = currentLevel;
		node_list.push_back(portalNode);

		while (1)
		{
			if ( node_list.count() == 0)
				break;

			currentLevel++;
			qDebug() << "node level " << currentLevel;

			QList<BranchNode_t*> tmp_list;
			tmp_list = node_list;
			node_list.clear();
			for (int i = 0; i < tmp_list.count(); i++)
			{
				for (int j = 0; j < tmp_list[i]->Children.count(); j++)
				{
					tmp_list[i]->Children[j]->depthFromPortal = currentLevel;
					node_list.push_back(tmp_list[i]->Children[j]);
				}
			}
		}
		

	}





	template<class TPixelType>
	void createVesselSeparation()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;

		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();
		typename TImageType::Pointer m_parentImage = ((VolumeSource<TImageType>*)currentSource()->parent())->getVolume();
		if (!m_parentImage)
		{
			cout << "Distance map is missing, cannot perform separation" << endl;
			return;
		}

		// Create a new volume with removed control points
		typename ImageDuplicator::Pointer dup0 = ImageDuplicator::New();
		dup0->SetInputImage(m_inputImage);
		dup0->Update();
		dup0->GetOutput()->FillBuffer(0);
		typename TImageType::Pointer m_outputImage = dup0->GetOutput();

		typename ImageDuplicator::Pointer dup1 = ImageDuplicator::New();
		dup1->SetInputImage(m_inputImage);
		dup1->Update();
		dup1->GetOutput()->FillBuffer(0);
		typename TImageType::Pointer m_portalVeinVolume = dup1->GetOutput();

		typename ImageDuplicator::Pointer dup2 = ImageDuplicator::New();
		dup2->SetInputImage(m_inputImage);
		dup2->Update();
		dup2->GetOutput()->FillBuffer(0);
		typename TImageType::Pointer m_hepaticVeinVolume = dup2->GetOutput();

		typename TImageType::SizeType inputImageSize = m_inputImage->GetBufferedRegion().GetSize();

		typename TImageType::SpacingType spacing = m_inputImage->GetSpacing();
		cout << "spacing " << spacing[0] << ", " << spacing[1] << ", " << spacing[2] << endl;

		QVector4D startPointV4;

		//QStringList seedlist = m_pickingCallBackForVesselSeparation->getPickedPosition().split(" ");
		QStringList seedlist = m_ui.m_cbSeedList->currentText().split(" ");
		if ( seedlist.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return;
		}

		int startingSeed[3];
		double startingRadius = seedlist[3].toDouble();// (seedlist[3].toDouble()-0.5f)*5.0f;

		startingSeed[0] = seedlist[0].toUInt();//int(seedlist[0].toDouble()/spacing[0]);
		startingSeed[1] = seedlist[1].toUInt();//int(seedlist[1].toDouble()/spacing[1]);
		startingSeed[2] = seedlist[2].toUInt();//int(seedlist[2].toDouble()/spacing[2]);

		startingRadius = evaluateBranchRadius<TPixelType>(startingSeed, startingRadius);
		if ( startingRadius == 0)
		{
			printf("radius evaluation failed\n");
		}

		startPointV4.setX(int(seedlist[0].toUInt()));
		startPointV4.setY(int(seedlist[1].toUInt()));
		startPointV4.setZ(int(seedlist[2].toUInt()));
		startPointV4.setW(startingRadius);//;

		cout << "Picked point [" << startPointV4.x() << "," << startPointV4.y() << "," 
			<< startPointV4.z() << "," << startPointV4.w() << "]" << endl;

		QStack<QStack<QVector4D>*> mainStack;
		QStack<QVector4D> *firstStack = new QStack<QVector4D>;
		firstStack->push(startPointV4);
		firstStack->push(startPointV4);
		mainStack.push(firstStack);

		QStack<QVector4D> withdrawStack;

		while(1)
		{
			if ( mainStack.count() == 0)	// no control point left
				break;

			QStack<QVector4D>* currentControlStack = mainStack.last();	// get the last one

			if ( currentControlStack->count()<=1)	// if no child left, then pop it, continue
			{
				currentControlStack = mainStack.pop();
				delete currentControlStack;

				continue;
			}

			typename TImageType::IndexType centerPointIndex;
			QVector4D v4 = currentControlStack->pop();		// pop the last point of current stack

			centerPointIndex[0] = v4.x();
			centerPointIndex[1] = v4.y();
			centerPointIndex[2] = v4.z();

			int adjacentNum = 0;			// get number of adjacent points
			bool hasMoreBranches = false;	// it is a new control point or not

			QStack<QVector4D> *newControlPointStack = new QStack<QVector4D>;	// assume it is a new control point
			newControlPointStack->push(v4);

			for ( unsigned int l = centerPointIndex[0]-1; l <= centerPointIndex[0]+1; l++)
			{
				for ( unsigned int m = centerPointIndex[1]-1; m <= centerPointIndex[1]+1; m++ )
				{
					for ( unsigned int n = centerPointIndex[2]-1; n <= centerPointIndex[2]+1; n++)
					{
						if (l<=5||m<=5||n<=5||l>=inputImageSize[0]-5||m>=inputImageSize[1]-5||n>=inputImageSize[2]-5)
							continue;

						typename TImageType::IndexType neighbourPointIndex;
						neighbourPointIndex[0] = l;
						neighbourPointIndex[1] = m;
						neighbourPointIndex[2] = n;

						if ( m_inputImage->GetPixel(neighbourPointIndex) == 1)	// has child
						{
							if ( m_outputImage->GetPixel(neighbourPointIndex) == 0)	// not processed
							{
								double radius = (double)m_pDistanceMap->GetPixel(neighbourPointIndex);

								if ( l==centerPointIndex[0] && m == centerPointIndex[1] && n == centerPointIndex[2])
									;
								else
									newControlPointStack->push(QVector4D(l,m,n, radius));	// add to child of a control point

								hasMoreBranches = true;		// we say for this possible control
							}

							adjacentNum++;
						}
					}
				}
			}

			m_outputImage->SetPixel(centerPointIndex, 1);

			if ( adjacentNum >= 4)		// is a new control point
			{
				static int controlCount = 1;
				//printf("control count %d\n", controlCount);
				controlCount++;


				// check whether it is control point
				bool isControlPoint = false;
				std::vector<Node>::iterator cpIt = m_vectorControlPoints.begin();
				for (; cpIt != m_vectorControlPoints.end(); cpIt++)
				{
					if ( cpIt->x == v4.x() && cpIt->y == v4.y() && cpIt->z == v4.z())
					{
						isControlPoint = true;
					}
				}


				QVector4D currentControlPoint = currentControlStack->front();
				QVector4D newControlPoint = newControlPointStack->front();

				if ( isControlPoint == false || 
					(std::abs(currentControlPoint.x() - newControlPoint.x()) <= 2 &&
					std::abs(currentControlPoint.y() - newControlPoint.y()) <= 2 &&
					std::abs(currentControlPoint.z() - newControlPoint.z()) <= 2))
				{
					QStack<QVector4D>::Iterator it = newControlPointStack->begin();
					for ( ++it; it != newControlPointStack->end(); it++)
					{
						currentControlStack->push(*it);
					}

					delete newControlPointStack;
					continue;
				}

				int currentCP[3];
				currentCP[0] = currentControlPoint.x();
				currentCP[1] = currentControlPoint.y();
				currentCP[2] = currentControlPoint.z();
				//double radius = evaluateBranchRadius<TPixelType>(currentCP, currentControlPoint.w());


				if ( currentControlPoint.w()+0.2 > newControlPoint.w()) // this control point has branches
				{
					mainStack.push(newControlPointStack);
				}
				else
				{
					;
				}
			}
			else // not a new control point, then move points to old control point
			{
				QStack<QVector4D>::Iterator it = newControlPointStack->begin();
				for ( ++it; it != newControlPointStack->end(); it++)
				{
					currentControlStack->push(*it);
				}

				delete newControlPointStack;

				m_outputImage->SetPixel(centerPointIndex, 1);
			}
		}


		// get hepatic

		typename TImageType::IndexType tmpIndex;
		typename TImageType::IndexType tmpIndex2;

		for(unsigned int i = 0; i < inputImageSize[0]; i++)
		{
			for (unsigned int j = 0; j < inputImageSize[1]; j++)
			{
				for (unsigned int k = 0; k < inputImageSize[2]; k++)
				{
					tmpIndex[0] = i;
					tmpIndex[1] = j;
					tmpIndex[2] = k;


					// get portal vein

					if ( m_outputImage->GetPixel(tmpIndex) == 1)
					{
						float radius = m_pDistanceMap->GetPixel(tmpIndex);
						radius += 3;

						int l = tmpIndex[0];
						int m = tmpIndex[1];
						int n = tmpIndex[2];

						for(unsigned int i = l-radius; i < l+radius; i++)
						{
							for (unsigned int j = m - radius; j < m+radius; j++)
							{
								for (unsigned int k = n - radius; k < n+radius; k++)
								{
									if (i<=0||j<=0||k<=0||i>inputImageSize[0]-1||j>inputImageSize[1]-1||k>=inputImageSize[2]-1)
										continue;

									if ( (i-l)*(i-l) + (j-m)*(j-m) + (k-n)*(k-n) > radius*radius)
										continue;

									tmpIndex2[0] = i;
									tmpIndex2[1] = j;
									tmpIndex2[2] = k;

									if ( m_parentImage->GetPixel(tmpIndex2) !=0 )
										m_portalVeinVolume->SetPixel(tmpIndex2, 1);
								}
							}
						}
					}
					else if ( m_outputImage->GetPixel(tmpIndex) == 2)
						m_outputImage->SetPixel(tmpIndex, 0);

					if ( m_parentImage->GetPixel(tmpIndex) !=0 && m_portalVeinVolume->GetPixel(tmpIndex) == 0)
						m_hepaticVeinVolume->SetPixel(tmpIndex, 1);

				}
			}
		}

		if ( 1)
		{
			for(unsigned int i = 0; i < inputImageSize[0]; i++)
			{
				for (unsigned int j = 0; j < inputImageSize[1]; j++)
				{
					for (unsigned int k = 0; k < inputImageSize[2]; k++)
					{
						tmpIndex[0] = i;
						tmpIndex[1] = j;
						tmpIndex[2] = k;

						if ( m_hepaticVeinVolume->GetPixel(tmpIndex) ==1)
						{
							int areaCount = 0;
							QVector4D v4;
							QStack<QVector4D> pointStack;

							v4.setX(i);	v4.setY(j);	v4.setZ(k);	v4.setW(0);
							pointStack.push(v4);

							while(1)
							{
								if ( pointStack.count()==0)
									break;

								typename TImageType::IndexType tmpIndex;
								QVector4D topPoint = pointStack.pop();
								tmpIndex[0] = topPoint.x();
								tmpIndex[1] = topPoint.y();
								tmpIndex[2] = topPoint.z();

								for ( unsigned int l = tmpIndex[0]-1; l <= tmpIndex[0]+1; l++)
								{
									for ( unsigned int m = tmpIndex[1]-1; m <= tmpIndex[1]+1; m++ )
									{
										for ( unsigned int n = tmpIndex[2]-1; n <= tmpIndex[2]+1; n++)
										{
											typename TImageType::IndexType tmpIndex2;

											tmpIndex2[0] = l;
											tmpIndex2[1] = m;
											tmpIndex2[2] = n;

											if (l<=1||m<=1||n<=1||l>=inputImageSize[0]-2||m>=inputImageSize[1]-2||n>=inputImageSize[2]-2)
												continue;

											if ( m_hepaticVeinVolume->GetPixel(tmpIndex2) == 1)
											{
												areaCount++;
												pointStack.push(QVector4D(l,m,n, 0));
												m_hepaticVeinVolume->SetPixel(tmpIndex2, 200); //temp 200
											}
										}
									}
								}
							}

							if ( areaCount <= 5000)
							{
								pointStack.clear();
								pointStack.push(v4);

								while(1)
								{
									if ( pointStack.count()==0)
										break;

									typename TImageType::IndexType tmpIndex;
									QVector4D topPoint = pointStack.pop();
									tmpIndex[0] = topPoint.x();
									tmpIndex[1] = topPoint.y();
									tmpIndex[2] = topPoint.z();

									for ( unsigned int l = tmpIndex[0]-1; l <= tmpIndex[0]+1; l++)
									{
										for ( unsigned int m = tmpIndex[1]-1; m <= tmpIndex[1]+1; m++ )
										{
											for ( unsigned int n = tmpIndex[2]-1; n <= tmpIndex[2]+1; n++)
											{
												typename TImageType::IndexType tmpIndex2;

												tmpIndex2[0] = l;
												tmpIndex2[1] = m;
												tmpIndex2[2] = n;

												if (l<=0||m<=0||n<=0||l>=inputImageSize[0]-1||m>=inputImageSize[1]-1||n>=inputImageSize[2]-1)
													continue;

												if ( m_hepaticVeinVolume->GetPixel(tmpIndex2) == 200)
												{
													pointStack.push(QVector4D(l,m,n, 0));
													m_hepaticVeinVolume->SetPixel(tmpIndex2, 0); //temp 200
												}
											}
										}
									}
								}
							}
						}
					}
				}
			}
		}



		VolumeSource<TImageType> *subsource = 
			new VolumeSource<TImageType>(QString("PortalSkeleton#%1").arg(currentSource()->rowCount()+1));
		subsource->setVolume(m_outputImage);
		currentSource()->addChild(subsource);


		VolumeSource<TImageType> *subsource1 = 
			new VolumeSource<TImageType>(QString("PortalVein#%1").arg(currentSource()->rowCount()+1));
		subsource1->setVolume(m_portalVeinVolume);
		currentSource()->addChild(subsource1);

		VolumeSource<TImageType> *subsource2 = 
			new VolumeSource<TImageType>(QString("HepaticVein#%1").arg(currentSource()->rowCount()+1));
		subsource2->setVolume(m_hepaticVeinVolume);
		currentSource()->addChild(subsource2);
	}


	

	template<class TPixelType>
	Source* createSmoothVolume(Source* source, int intenal)
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::Image<unsigned char, 3> TImageUChar3D;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::Image<float, 3> RealImageType;

		typedef itk::CastImageFilter<TImageType, RealImageType> CastToRealFilterType;
		typedef itk::AntiAliasBinaryImageFilter<RealImageType, RealImageType> FilterType;
		typedef itk::RescaleIntensityImageFilter<RealImageType, TImageUChar3D> RescaleFilter;
		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)source)->getVolume();

		typename CastToRealFilterType::Pointer toReal = CastToRealFilterType::New();
		typename RescaleFilter::Pointer rescale = RescaleFilter::New();
		typename FilterType::Pointer filter = FilterType::New();

		toReal->SetInput(m_inputImage);
		toReal->Update();

		filter->SetInput(toReal->GetOutput());
		if (intenal)
			filter->SetNumberOfIterations(10);
		else
			filter->SetNumberOfIterations(m_ui.m_textSMVolumeIteration->text().toInt());
		filter->SetMaximumRMSError(0.01);
		filter->SetNumberOfLayers(2);
		filter->Update();

		rescale->SetInput(filter->GetOutput());
		rescale->Update();

		typename TImageUChar3D::IndexType indexType;
		typename TImageUChar3D::SizeType size = rescale->GetOutput()->GetBufferedRegion().GetSize();
		for (unsigned int i = 0; i < size[0]; i++)
		{
			for ( unsigned int j = 0; j < size[1]; j++)
			{
				for ( unsigned int k = 0; k < size[2]; k++)
				{
					indexType[0] = i;
					indexType[1] = j;
					indexType[2] = k;

					if ( rescale->GetOutput()->GetPixel(indexType) > 0)
						rescale->GetOutput()->SetPixel(indexType, 1);
				}
			}
		}


		
		VolumeSource<TImageUChar3D> *subsource = 
			new VolumeSource<TImageUChar3D>(QString("SmoothVolume#%1").arg(currentSource()->rowCount()+1));
		subsource->setVolume(rescale->GetOutput());


// 		VolumeSource<RealImageType> *subsource = 
// 			new VolumeSource<RealImageType>(QString("SmoothVolume#%1").arg(currentSource()->rowCount()+1));
// 		subsource->setVolume(filter->GetOutput());

		if ( intenal == 0)
			currentSource()->addChild(subsource);

		return (Source*)subsource;
	}


	template<class TImageType> 
	void divideLiverBySeedWithDirections( int add2ui, QList<QVector4D> seeds, QList<QVector4D> seedDirections, 
						typename TImageType::Pointer volumeSkeleton, typename TImageType::Pointer volumeLiver, int color )
	{
		typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::RegionOfInterestImageFilter<TImageType, TImageType> RegionOfInterestImageFilter;
		typedef itk::SignedDanielssonDistanceMapImageFilter<TImageType, TImageType> SignedDanielssonDistanceMapImageFilter;

		itk::TimeProbe clock; clock.Start();
		printf("[time] dividing liver by selected branches...\n");

		typename ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
		maskSO->SetImage(volumeLiver);
		maskSO->Update();

		typename TImageType::RegionType regionType = maskSO->GetAxisAlignedBoundingBoxRegion();
// 		typename TImageType::IndexType regionIndex = regionType.GetIndex();
// 		typename TImageType::SizeType regionSize = regionType.GetSize();
// 
// 		printf("original region info (%d,%d,%d) (%d,%d,%d)\n", regionIndex[0], regionIndex[1], regionIndex[2], regionSize[0], regionSize[1], regionSize[2]);
// 		regionIndex[0] -= 30;
// 		regionIndex[1] -= 30;
// 		regionIndex[2] -= 30;
// 		regionSize[0] += 60;
// 		regionSize[1] += 60;
// 		regionSize[2] += 60;
// 		regionType.SetIndex(regionIndex);
// 		regionType.SetSize(regionSize);
// 
// 		printf("current region info (%d,%d,%d) (%d,%d,%d)\n", regionIndex[0], regionIndex[1], regionIndex[2], regionSize[0], regionSize[1], regionSize[2]);

		typename RegionOfInterestImageFilter::Pointer roiLiverFilter = RegionOfInterestImageFilter::New();
		roiLiverFilter->SetInput(volumeLiver);
		roiLiverFilter->SetRegionOfInterest(regionType);
		roiLiverFilter->Update();

		typename RegionOfInterestImageFilter::Pointer roiSKFilter = RegionOfInterestImageFilter::New();
		maskSO->SetImage(volumeSkeleton);
		maskSO->Update();
		roiSKFilter->SetInput(volumeSkeleton);
		roiSKFilter->SetRegionOfInterest(regionType);
		roiSKFilter->Update();

		typename TImageType::Pointer roiLiver = roiLiverFilter->GetOutput();
		typename TImageType::Pointer roiSkeleton = roiSKFilter->GetOutput();

		typename ImageDuplicator::Pointer dupForSelectedBranches = ImageDuplicator::New();
		dupForSelectedBranches->SetInputImage(roiSkeleton); 
		dupForSelectedBranches->Update();
		typename ImageDuplicator::Pointer dupForUnSelectedBranches = ImageDuplicator::New();
		dupForUnSelectedBranches->SetInputImage(roiSkeleton); 
		dupForUnSelectedBranches->Update();

		typename TImageType::Pointer pointSelectedBranches = dupForSelectedBranches->GetOutput();
		typename TImageType::Pointer pointUnSelectedBranches = dupForUnSelectedBranches->GetOutput();
		pointSelectedBranches->FillBuffer(0);

		typename TImageType::IndexType seedIndex;
		typename TImageType::IndexType directionIndex;

		for (int i = 0; i < seeds.count(); i++)
		{
			seedIndex[0] = seeds.at(i).x() - regionType.GetIndex()[0];
			seedIndex[1] = seeds.at(i).y() - regionType.GetIndex()[1];
			seedIndex[2] = seeds.at(i).z() - regionType.GetIndex()[2];

			dupForSelectedBranches->GetOutput()->SetPixel(seedIndex, 1);

			directionIndex[0] = seedDirections.at(i).x() - regionType.GetIndex()[0];
			directionIndex[1] = seedDirections.at(i).y() - regionType.GetIndex()[1];
			directionIndex[2] = seedDirections.at(i).z() - regionType.GetIndex()[2];

			QStack<typename TImageType::IndexType> indexStack;
			indexStack.push(directionIndex);

			while(1)
			{
				if ( indexStack.size() == 0)
					break;

				typename TImageType::IndexType currIndex;
				currIndex = indexStack.pop();

				for ( int i = currIndex[0]-1; i <= currIndex[0]+1; i++)
				{
					for ( int j = currIndex[1]-1; j <= currIndex[1]+1; j++ )
					{
						for ( int k = currIndex[2]-1; k <= currIndex[2]+1; k++)
						{
							typename TImageType::IndexType tmpIndex;

							tmpIndex[0] = i;
							tmpIndex[1] = j;
							tmpIndex[2] = k;

							if ( roiSkeleton->GetPixel(tmpIndex) != 0)		// in branch
							{
								if ( pointSelectedBranches->GetPixel(tmpIndex) == 0)	//
								{
									indexStack.push(tmpIndex);
									pointSelectedBranches->SetPixel(tmpIndex, 1);
									pointUnSelectedBranches->SetPixel(tmpIndex, 0);
								}
							}
						}
					}
				}
			}
		}



		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime()); clock.Start();


		typedef itk::DanielssonDistanceMapImageFilter<TImageType, TImageType> DanielssonDistanceMapImageFilter;

		typename DanielssonDistanceMapImageFilter::Pointer dmFilter1 = DanielssonDistanceMapImageFilter::New();
		dmFilter1->SetInput(pointSelectedBranches);
		//dmFilter1->InputIsBinaryOn();
		dmFilter1->Update();

		typename DanielssonDistanceMapImageFilter::Pointer dmFilter2 = DanielssonDistanceMapImageFilter::New();
		dmFilter2->SetInput(pointUnSelectedBranches);
		//dmFilter2->InputIsBinaryOn();
		dmFilter2->Update();

// 		VolumeSource<TImageType> *subsource1 = 
// 			new VolumeSource<TImageType>(QString("Selected"));
// 		subsource1->setVolume(dmFilter1->GetOutput());
// 		m_sourceLiver->addChild(subsource1);
// 
// 		VolumeSource<TImageType> *subsource2 = 
// 			new VolumeSource<TImageType>(QString("UnSelected"));
// 		subsource2->setVolume(dmFilter2->GetOutput());
// 		m_sourceLiver->addChild(subsource2);

		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime()); clock.Start();
		printf("[time] separating liver volume...\n");

		typename ImageDuplicator::Pointer dupLiver1 = ImageDuplicator::New();
		dupLiver1->SetInputImage(volumeSkeleton); 
		dupLiver1->Update();
		dupLiver1->GetOutput()->FillBuffer(0);

		typename ImageDuplicator::Pointer dupLiver2 = ImageDuplicator::New();
		dupLiver2->SetInputImage(dupLiver1->GetOutput()); 
		dupLiver2->Update();

		typename ImageDuplicator::Pointer dupCuttingPanel = ImageDuplicator::New();
		dupCuttingPanel->SetInputImage(dupLiver1->GetOutput()); 
		dupCuttingPanel->Update();

		typename TImageType::SizeType size = pointSelectedBranches->GetBufferedRegion().GetSize();
		typename TImageType::IndexType indexType;
		typename TImageType::IndexType origIndexType;
		typename TImageType::Pointer dmPointer1 = dmFilter1->GetOutput();
		typename TImageType::Pointer dmPointer2 = dmFilter2->GetOutput();
		typename TImageType::Pointer dupPointer1 = dupLiver1->GetOutput();
		typename TImageType::Pointer dupPointer2 = dupLiver2->GetOutput();
		typename TImageType::Pointer dupCuttingPanelPoint = dupCuttingPanel->GetOutput();
		typename TImageType::PixelType pixel1;
		typename TImageType::PixelType pixel2;

		vtkPoints *selectionPoints = vtkPoints::New();

		for ( unsigned int i = 0; i < size[0]; i++ )
		{
			for ( unsigned int j = 0; j < size[1]; j++)
			{
				for ( unsigned int k = 0; k < size[2]; k++)
				{
					indexType[0] = i;
					indexType[1] = j;
					indexType[2] = k;

					if ( roiLiver->GetPixel(indexType) == 0)
						continue;

					origIndexType[0] = indexType[0]+regionType.GetIndex()[0];
					origIndexType[1] = indexType[1]+regionType.GetIndex()[1];
					origIndexType[2] = indexType[2]+regionType.GetIndex()[2];

					pixel1 = dmPointer1->GetPixel(indexType);
					pixel2 = dmPointer2->GetPixel(indexType);

					if( pixel1 < pixel2)
					{
						dupPointer1->SetPixel(origIndexType, 1);

					}
					else 
					{
						dupPointer2->SetPixel(origIndexType, 1);
					}


					if( pixel1 == pixel2)
					{
// 						selectionPoints->InsertNextPoint((double)origIndexType[0]/m_currentSpacing[0], 
// 														 (double)origIndexType[1]/m_currentSpacing[1], 
// 														 (double)origIndexType[2]/m_currentSpacing[2]);
						dupCuttingPanelPoint->SetPixel(origIndexType, 1);
					}
				}
			}
		}




		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime());

		if ( add2ui)
		{
			VolumeSource<TImageType> *subsource1 = 
				new VolumeSource<TImageType>(QString("LiverPart"));
			subsource1->setVolume(dupLiver1->GetOutput());

			Source* s = createSmoothVolume<TImageType::PixelType>(subsource1, 1);
			SurfaceSource* s2 = volume2SmoothSurface(s, QString("Liverpart"), 1);
			s2->sourceModel();
			s2->setSurfaceColor(color);

			VolumeSource<TImageType> *cuttingPanel = 
				new VolumeSource<TImageType>(QString("CuttingPanel"));
			cuttingPanel->setVolume(dupCuttingPanel->GetOutput());

			m_sourceLiver->addChild(s2);
			m_sourceLiver->addChild(subsource1);
			m_sourceLiver->addChild(cuttingPanel);


// 			SurfaceSource* liverSurface = volume2SmoothSurface(m_sourceLiver, QString("LiverSurface"), 1);
// 			liverSurface->sourceModel();
// 			m_sourceLiver->addChild(liverSurface);
// 
// 			vtkSelectPolyData *loop = vtkSelectPolyData::New();
// 			loop->SetInput(liverSurface->sdata());
// 			loop->SetLoop(selectionPoints);
// 			loop->GenerateSelectionScalarsOn();
// 			loop->SetSelectionModeToSmallestRegion();
// 
// 			vtkClipPolyData *clip = vtkClipPolyData::New();
// 			clip->SetInputConnection(loop->GetOutputPort());
// 
// 			vtkPolyDataMapper *clipMapper = vtkPolyDataMapper::New();
// 			clipMapper->SetInputConnection(clip->GetOutputPort());
// 
// 			vtkLODActor *clipActor = vtkLODActor::New();
// 			clipActor->SetMapper(clipMapper);
// 
// 			m_displayModel->addActor(clipActor);

		}
	}


	template<class TPixelType>
	void testSmoothing()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typename TImageType::Pointer m_inputImage = ((VolumeSource<TImageType>*)currentSource())->getVolume();

		VolumeSource<TImageType> *subsource1 = 
			new VolumeSource<TImageType>(QString("LiverPart"));
		subsource1->setVolume(m_inputImage);


 		Source* s = createSmoothVolume<TImageType::PixelType>(subsource1, 1);
 		SurfaceSource* s2 = volume2SmoothSurface(s, QString("Liverpart"), 1);

		//s2->setSurfaceColor(3);

		currentSource()->addChild(s2);

	}




	template<class TImageType>
	void divideLiverBySeed(int add2ui, typename TImageType::IndexType seed, 
							typename TImageType::IndexType seedDirection,
							typename TImageType::Pointer volumeSkeleton,
							typename TImageType::Pointer volumeLiver,
							typename TImageType::Pointer *outputSkeleton1,
							typename TImageType::Pointer *outputSkeleton2,
							typename TImageType::Pointer *outputLiver1,
							typename TImageType::Pointer *outputLiver2,
							typename TImageType::IndexType *seed1,
							typename TImageType::IndexType *seed2)
	{
 		typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;
		typedef itk::ImageDuplicator<TImageType> ImageDuplicator;
		typedef itk::RegionOfInterestImageFilter<TImageType, TImageType> RegionOfInterestImageFilter;
		typedef itk::SignedDanielssonDistanceMapImageFilter<TImageType, TImageType> SignedDanielssonDistanceMapImageFilter;

		itk::TimeProbe clock; clock.Start();
		printf("[time] separating vessels...\n");

		typename ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
		maskSO->SetImage(volumeLiver);
		maskSO->Update();

		typename TImageType::RegionType regionType = maskSO->GetAxisAlignedBoundingBoxRegion();

		typename RegionOfInterestImageFilter::Pointer roiLiverFilter = RegionOfInterestImageFilter::New();
		roiLiverFilter->SetInput(volumeLiver);
		roiLiverFilter->SetRegionOfInterest(regionType);
		roiLiverFilter->Update();

		typename RegionOfInterestImageFilter::Pointer roiSKFilter = RegionOfInterestImageFilter::New();
		maskSO->SetImage(volumeSkeleton);
		maskSO->Update();
		roiSKFilter->SetInput(volumeSkeleton);
		roiSKFilter->SetRegionOfInterest(regionType);
		roiSKFilter->Update();

		typename TImageType::Pointer roiLiver = roiLiverFilter->GetOutput();
		typename TImageType::Pointer roiSkeleton = roiSKFilter->GetOutput();
		typename TImageType::IndexType tmpIndex;
		typename TImageType::IndexType initialIndex;

		initialIndex[0] = seed[0] - regionType.GetIndex()[0];
		initialIndex[1] = seed[1] - regionType.GetIndex()[1];
		initialIndex[2] = seed[2] - regionType.GetIndex()[2];

		QList<typename TImageType::IndexType> branchIndexList;
		QList<typename TImageType::IndexType> branchNewSeedList;
		QList<typename TImageType::Pointer>   branchVolumeList;
		

		int showNecessaryDivision = 0;
		cout << "---- separate branches based on seed point [ "
			 << seed[0] << ", "
			 << seed[1] << ", "
			 << seed[2] << "]" << endl;
		for ( int i = initialIndex[0]-1; i <= initialIndex[0]+1; i++)
		{
			for ( int j = initialIndex[1]-1; j <= initialIndex[1]+1; j++ )
			{
				for ( int k = initialIndex[2]-1; k <= initialIndex[2]+1; k++)
				{
					tmpIndex[0] = i;
					tmpIndex[1] = j;
					tmpIndex[2] = k;

					if ( i == initialIndex[0] && j == initialIndex[1] && k == initialIndex[2])
						continue;

					if ( roiSkeleton->GetPixel(tmpIndex) != 0)
					{
						static int branch = 1;
						cout << "  create branch " << branch++ << " :[" << i + regionType.GetIndex()[0] << ", " 
							                                     << j + regionType.GetIndex()[1] << ", " 
																 << k + regionType.GetIndex()[2] << "]" << endl;

						branchIndexList.push_back(tmpIndex);

						typename ImageDuplicator::Pointer dup = ImageDuplicator::New();
						dup->SetInputImage(roiSkeleton); 
						dup->Update();
						dup->GetOutput()->FillBuffer(0);
						dup->GetOutput()->SetPixel(tmpIndex, 1);
						dup->GetOutput()->SetPixel(initialIndex, 1);
						branchVolumeList.push_back(dup->GetOutput());

						//tmpIndex[0] = 0; tmpIndex[1]=0; tmpIndex[2]=0;
						branchNewSeedList.push_back(tmpIndex);
					}
				}
			}
		}

		// fill the square which contains starting point to stop growing back
		for ( int i = initialIndex[0]-1; i <= initialIndex[0]+1; i++)
		{
			for ( int j = initialIndex[1]-1; j <= initialIndex[1]+1; j++ )
			{
				for ( int k = initialIndex[2]-1; k <= initialIndex[2]+1; k++)
				{
					tmpIndex[0] = i;
					tmpIndex[1] = j;
					tmpIndex[2] = k;

					typename QList<typename TImageType::Pointer>::Iterator it = branchVolumeList.begin();
					for (;it != branchVolumeList.end(); it++)
					{
						(*it)->SetPixel(tmpIndex, 1);
					}
				}
			}
		}


		cout << "---- find out seed point (the next point of starting point) of each branch ----" << endl;
		QStack<typename TImageType::IndexType> indexStack;
		QList<QList<typename TImageType::IndexType>*> indexLists;

		for (int branch = 0; branch < branchIndexList.size(); branch++)
		{
			indexStack.clear();
			indexStack.push(branchIndexList[branch]);
            QList<typename TImageType::IndexType>* indexlist = new QList<typename TImageType::IndexType>();
			indexLists.push_back(indexlist);

			typename TImageType::IndexType nextIndex;
			typename TImageType::Pointer dupVolume = branchVolumeList[branch];
			int gotTheFirstSeed = 0;

			while (1)
			{
				if (indexStack.size() == 0)
					break;

				nextIndex =	indexStack.pop();

				int radius = 1;
				int neighborCount = 0;
				for ( int i = nextIndex[0]-radius; i <= nextIndex[0]+radius; i++)
				{
					for ( int j = nextIndex[1]-radius; j <= nextIndex[1]+radius; j++ )
					{
						for ( int k = nextIndex[2]-radius; k <= nextIndex[2]+radius; k++)
						{
							tmpIndex[0] = i;
							tmpIndex[1] = j;
							tmpIndex[2] = k;

							if ( roiSkeleton->GetPixel(tmpIndex) != 0)
							{
								neighborCount++;

								if (dupVolume->GetPixel(tmpIndex)==0)
								{
									dupVolume->SetPixel(tmpIndex, 1);

									indexStack.push(tmpIndex);
									indexlist->push_back(tmpIndex);
								}
							}
						}
					}
				}

				if ( neighborCount >= 4 && gotTheFirstSeed == 0)
				{
					if (std::labs(branchNewSeedList[branch][0] - nextIndex[0]) >= 3 || 
						std::labs(branchNewSeedList[branch][1] - nextIndex[1]) >= 3 || 
						std::labs(branchNewSeedList[branch][2] - nextIndex[2]) >= 3)
					{
						branchNewSeedList[branch][0] = nextIndex[0];
						branchNewSeedList[branch][1] = nextIndex[1];
						branchNewSeedList[branch][2] = nextIndex[2];
						gotTheFirstSeed = 1;

						cout << "  got next seed for branch " << branch << " : "
							<< nextIndex[0] + regionType.GetIndex()[0] << ", " 
							<< nextIndex[1] + regionType.GetIndex()[1] << ", " 
							<< nextIndex[2] + regionType.GetIndex()[2] << endl;
					}
				}
			}
// 
// 
// 
// 			VolumeSource<TImageType> *subsource1 = 
// 				new VolumeSource<TImageType>(QString("Branches"));
// 			subsource1->setVolume(dupVolume);
// 			m_sourceSkeleton->addChild(subsource1);
		}

		
		// remove duplicated seed list
		for (int i = 0; i < indexLists.size(); i++)
		{
			for (int j = i+1; j < indexLists.size(); j++)
			{
				if ( indexLists[i]->size() == indexLists[j]->size())
				{
					indexLists[j]->clear();
				}
			}
		}



		cout << "---- remove duplicated index list ----" << endl;
		while(indexLists.size()>2)
		{
			int smallestIndex = 0;
			int smallestSize = 100000;

			typename QList<QList<typename TImageType::IndexType>*>::Iterator it = indexLists.begin();
			for (int i = 0; it != indexLists.end(); i++, it++)
			{
				if ( (*it)->size() < smallestSize )
				{
					smallestIndex = i;
					smallestSize = (*it)->size();
				}
			}

			cout << "  remove index list " << smallestIndex << " which has " << smallestSize << " points" << endl;
			cout << "  remove seed point [" << branchNewSeedList[smallestIndex][0] + regionType.GetIndex()[0] << ", " 
										    << branchNewSeedList[smallestIndex][1] + regionType.GetIndex()[1] << ", "
										    << branchNewSeedList[smallestIndex][2] + regionType.GetIndex()[2] << "]" << endl;
			branchVolumeList.removeAt(smallestIndex);
			branchNewSeedList.removeAt(smallestIndex);
			indexLists.removeAt(smallestIndex);
		}

		cout << "  left index 1: size[" << indexLists[0]->size() << "], next seed[" << branchNewSeedList[0][0] + regionType.GetIndex()[0]<< ","
																					<< branchNewSeedList[0][1] + regionType.GetIndex()[1]<< ","
																					<< branchNewSeedList[0][2] + regionType.GetIndex()[2]<< "]" << endl;
		cout << "  left index 2: size[" << indexLists[1]->size() << "], next seed[" << branchNewSeedList[1][0] + regionType.GetIndex()[0]<< ","
																					<< branchNewSeedList[1][1] + regionType.GetIndex()[1]<< ","
																					<< branchNewSeedList[1][2] + regionType.GetIndex()[2]<< "]" << endl;

		if ( seed1 != NULL)
		{
			(*seed1)[0] = branchNewSeedList[0][0] + regionType.GetIndex()[0];
			(*seed1)[1] = branchNewSeedList[0][1] + regionType.GetIndex()[1];
			(*seed1)[2] = branchNewSeedList[0][2] + regionType.GetIndex()[2];

			(*seed2)[0] = branchNewSeedList[1][0] + regionType.GetIndex()[0];
			(*seed2)[1] = branchNewSeedList[1][1] + regionType.GetIndex()[1];
			(*seed2)[2] = branchNewSeedList[1][2] + regionType.GetIndex()[2];
		}


		if ( indexLists.size() != 2)
		{
			cout << "indexlist.size() != 2" << endl;
			return;
		}
		else
		{
			typename TImageType::SizeType size = branchVolumeList[0]->GetBufferedRegion().GetSize();
			typename TImageType::IndexType indexType;
			typename TImageType::IndexType origIndexType;

			cout << "branch0[0]:" << (branchIndexList[0])[0]  + regionType.GetIndex()[0]
						   << "," << (branchIndexList[0])[1]  + regionType.GetIndex()[1]
						   << "," << (branchIndexList[0])[2]  + regionType.GetIndex()[2] << endl;
			cout << "branch1[0]:" << (branchIndexList[1])[0]  + regionType.GetIndex()[0]
			               << "," << (branchIndexList[1])[1]  + regionType.GetIndex()[1]
						   << "," << (branchIndexList[1])[2]  + regionType.GetIndex()[2] << endl;

			if ( (branchIndexList[0])[0] + regionType.GetIndex()[0] == seedDirection[0] && 
				 (branchIndexList[0])[1] + regionType.GetIndex()[1] == seedDirection[1] && 
				 (branchIndexList[0])[2] + regionType.GetIndex()[2] == seedDirection[2]) 
			{
				showNecessaryDivision = 1;
			}
			else if ( (branchIndexList[1])[0] + regionType.GetIndex()[0] == seedDirection[0] && 
				      (branchIndexList[1])[1] + regionType.GetIndex()[1] == seedDirection[1] && 
					  (branchIndexList[1])[2] + regionType.GetIndex()[2] == seedDirection[2]) 
			{
				showNecessaryDivision = 2;
			}
			else
				showNecessaryDivision = 0;



			typename ImageDuplicator::Pointer dupSkeleton1 = ImageDuplicator::New();
			dupSkeleton1->SetInputImage(volumeSkeleton); 
			dupSkeleton1->Update();
			dupSkeleton1->GetOutput()->FillBuffer(0);
            
			typename ImageDuplicator::Pointer dupSkeleton2 = ImageDuplicator::New();
			dupSkeleton2->SetInputImage(volumeSkeleton); 
			dupSkeleton2->Update();
			dupSkeleton2->GetOutput()->FillBuffer(0);

			typename TImageType::Pointer pDupSkeleton1 = dupSkeleton1->GetOutput();
			typename TImageType::Pointer pDupSkeleton2 = dupSkeleton2->GetOutput();

			for ( unsigned int i = 0; i < size[0]; i++ )
			{
				for ( unsigned int j = 0; j < size[1]; j++)
				{
					for ( unsigned int k = 0; k < size[2]; k++)
					{
						indexType[0] = i;
						indexType[1] = j;
						indexType[2] = k;

						if ( roiSkeleton->GetPixel(indexType) == 0)
							continue;

						origIndexType[0] = indexType[0]+regionType.GetIndex()[0];
						origIndexType[1] = indexType[1]+regionType.GetIndex()[1];
						origIndexType[2] = indexType[2]+regionType.GetIndex()[2];

						if( branchVolumeList[0]->GetPixel(indexType))
						{
							pDupSkeleton1->SetPixel(origIndexType, 1);
						}
						else
						{
							pDupSkeleton2->SetPixel(origIndexType, 1);
						}
					}
				}
			}

// 			VolumeSource<TImageType> *subsource1 = 
// 				new VolumeSource<TImageType>(QString("Branch1"));
// 			subsource1->setVolume(pDupSkeleton1);
// 			m_sourceSkeleton->addChild(subsource1);
// 
// 			subsource1 = new VolumeSource<TImageType>(QString("Branch2"));
// 			subsource1->setVolume(pDupSkeleton2);
// 			m_sourceSkeleton->addChild(subsource1);

			if ( outputSkeleton1 != NULL)
			{
				*outputSkeleton1 = pDupSkeleton1;
				*outputSkeleton2 = pDupSkeleton2;
			}
		}

		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime()); clock.Start();
		printf("[time] generating distance map...\n");


		
// 		typedef itk::SignedMaurerDistanceMapImageFilter<TImageType, TImageType> SignedMaurerDistanceMapImageFilter;
// 		SignedMaurerDistanceMapImageFilter::Pointer dmFilter1 = SignedMaurerDistanceMapImageFilter::New();
// 		dmFilter1->SetInput(branchVolumeList[0]);
// 		dmFilter1->Update();
// 
// 		SignedMaurerDistanceMapImageFilter::Pointer dmFilter2 = SignedMaurerDistanceMapImageFilter::New();
// 		dmFilter2->SetInput(branchVolumeList[1]);
// 		dmFilter2->Update();

// 		typedef itk::ApproximateSignedDistanceMapImageFilter<TImageType, TImageType> ApproximateSignedDistanceMapImageFilter;
// 		ApproximateSignedDistanceMapImageFilter::Pointer dmFilter1 = ApproximateSignedDistanceMapImageFilter::New();
// 		dmFilter1->SetInput(branchVolumeList[0]);
// 		dmFilter1->Update();
// 
// 		ApproximateSignedDistanceMapImageFilter::Pointer dmFilter2 = ApproximateSignedDistanceMapImageFilter::New();
// 		dmFilter2->SetInput(branchVolumeList[1]);
// 		dmFilter2->Update();

		typedef itk::DanielssonDistanceMapImageFilter<TImageType, TImageType> DanielssonDistanceMapImageFilter;
		typename DanielssonDistanceMapImageFilter::Pointer dmFilter1 = DanielssonDistanceMapImageFilter::New();
		dmFilter1->SetInput(branchVolumeList[0]);
		dmFilter1->InputIsBinaryOn();
		dmFilter1->Update();

		typename DanielssonDistanceMapImageFilter::Pointer dmFilter2 = DanielssonDistanceMapImageFilter::New();
		dmFilter2->SetInput(branchVolumeList[1]);
		dmFilter2->InputIsBinaryOn();
		dmFilter2->Update();

// 		SignedDanielssonDistanceMapImageFilter::Pointer dmFilter1 = SignedDanielssonDistanceMapImageFilter::New();
// 		dmFilter1->SetInput(branchVolumeList[0]);
// 		dmFilter1->Update();
// 
// 		SignedDanielssonDistanceMapImageFilter::Pointer dmFilter2 = SignedDanielssonDistanceMapImageFilter::New();
// 		dmFilter2->SetInput(branchVolumeList[1]);
// 		dmFilter2->Update();

// 		VolumeSource<TImageType> *subsource3 = 
// 			new VolumeSource<TImageType>(QString("Distance1"));
// 		subsource3->setVolume(dmFilter1->GetOutput());
// 		m_sourceSkeleton->addChild(subsource3);
// 
// 		VolumeSource<TImageType> *subsource4 = 
// 			new VolumeSource<TImageType>(QString("Distance2"));
// 		subsource4->setVolume(dmFilter2->GetOutput());
// 		m_sourceSkeleton->addChild(subsource4);

		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime()); clock.Start();
		printf("[time] separating liver volume...\n");

		typename ImageDuplicator::Pointer dupLiver1 = ImageDuplicator::New();
		dupLiver1->SetInputImage(volumeSkeleton); 
		dupLiver1->Update();
		dupLiver1->GetOutput()->FillBuffer(0);

		typename ImageDuplicator::Pointer dupLiver2 = ImageDuplicator::New();
		dupLiver2->SetInputImage(dupLiver1->GetOutput()); 
		dupLiver2->Update();

		typename TImageType::SizeType size = branchVolumeList[0]->GetBufferedRegion().GetSize();
		typename TImageType::IndexType indexType;
		typename TImageType::IndexType origIndexType;

		for ( unsigned int i = 0; i < size[0]; i++ )
		{
			for ( unsigned int j = 0; j < size[1]; j++)
			{
				for ( unsigned int k = 0; k < size[2]; k++)
				{
					indexType[0] = i;
					indexType[1] = j;
					indexType[2] = k;

					if ( roiLiver->GetPixel(indexType) == 0)
						continue;

					origIndexType[0] = indexType[0]+regionType.GetIndex()[0];
					origIndexType[1] = indexType[1]+regionType.GetIndex()[1];
					origIndexType[2] = indexType[2]+regionType.GetIndex()[2];

					if( dmFilter1->GetOutput()->GetPixel(indexType) < dmFilter2->GetOutput()->GetPixel(indexType))
					{
						dupLiver1->GetOutput()->SetPixel(origIndexType, 1);
					}
					else
					{
						dupLiver2->GetOutput()->SetPixel(origIndexType, 1);
					}
				}
			}
		}
		

		clock.Stop(); printf("[done] : %.2fs\n", clock.GetMeanTime());

		if ( add2ui)
		{
			VolumeSource<TImageType> *subsource1 = 
				new VolumeSource<TImageType>(QString("LiverPart1"));
			subsource1->setVolume(dupLiver1->GetOutput());

			VolumeSource<TImageType> *subsource2 = 
				new VolumeSource<TImageType>(QString("LiverPart2"));
			subsource2->setVolume(dupLiver2->GetOutput());

			if ( showNecessaryDivision == 1)
			{
				m_sourceLiver->addChild(subsource1);
				m_sourceLiver = subsource2;
			}
			else if ( showNecessaryDivision == 2)
			{
				m_sourceLiver->addChild(subsource2);
				m_sourceLiver = subsource1;
			}
		}

		if ( outputLiver1 != NULL)
		{
			*outputLiver1 = dupLiver1->GetOutput();
			*outputLiver2 = dupLiver2->GetOutput();
		}
	}


	template<class TImageType>
	void evaluateSeed(typename TImageType::IndexType *seed)
	{
		cout << "---- evaluate seed [" << (*seed)[0] << "," 
									   << (*seed)[1] << ","
									   << (*seed)[2] << "] ----" << endl;
		for (int i = 0; i < m_ui.m_cbSeedList->count(); i++)
		{
			QStringList sl = (m_ui.m_cbSeedList->itemText(i)).split(" ");


			if (std::labs(sl[0].toUInt()-(int)((*seed)[0])) <= 4 && 
				std::labs(sl[1].toUInt()-(int)((*seed)[1])) <= 4 && 
				std::labs(sl[2].toUInt()-(int)((*seed)[2])) <= 4)
			{
				(*seed)[0] = sl[0].toUInt();
				(*seed)[1] = sl[1].toUInt();
				(*seed)[2] = sl[2].toUInt();
				cout << "  seed moved to : [" << (*seed)[0] << ", " << (*seed)[1] << ", " << (*seed)[2] << "]" << endl;
				break;
			}
		}
	}

	
	
	template<class TPixelType>
	void liverDivision()
	{
		typedef itk::Image<TPixelType, 3> TImageType;
		typedef itk::Image<unsigned char, 3> TImageUChar3D;
		typedef itk::CastImageFilter<TImageType, TImageUChar3D> CastImageFilter;
		typedef itk::ImageDuplicator<TImageUChar3D> ImageDuplicator;
		typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject;
		typedef itk::RegionOfInterestImageFilter<TImageUChar3D, TImageUChar3D> RegionOfInterestImageFilter;
		typedef itk::SignedDanielssonDistanceMapImageFilter<TImageUChar3D, TImageUChar3D> SignedDanielssonDistanceMapImageFilter;


		typename TImageType::Pointer m_volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		typename TImageType::Pointer m_volumeLiver = ((VolumeSource<TImageType>*)m_sourceLiver)->getVolume();


		QStringList seedlist1 = m_ui.m_cbSeedList->currentText().split(" ");
		if ( seedlist1.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return;
		}

		typename CastImageFilter::Pointer castFilter1 = CastImageFilter::New();
		castFilter1->SetInput(m_volumeSkeleton);
		castFilter1->Update();

		typename CastImageFilter::Pointer castFilter2 = CastImageFilter::New();
		castFilter2->SetInput(m_volumeLiver);
		castFilter2->Update();

		typename TImageType::SpacingType spacing = m_volumeLiver->GetSpacing();



		typename TImageUChar3D::IndexType seed;
		typename TImageUChar3D::Pointer outputSk11;
		typename TImageUChar3D::Pointer outputSk12;
		typename TImageUChar3D::Pointer outputSk21;
		typename TImageUChar3D::Pointer outputSk22;
		typename TImageUChar3D::Pointer outputSk23;
		typename TImageUChar3D::Pointer outputSk24;

		typename TImageUChar3D::Pointer outputLv11;
		typename TImageUChar3D::Pointer outputLv12;
		typename TImageUChar3D::Pointer outputLv21;
		typename TImageUChar3D::Pointer outputLv22;
		typename TImageUChar3D::Pointer outputLv23;
		typename TImageUChar3D::Pointer outputLv24;

		typename TImageUChar3D::IndexType outputSeed11;
		typename TImageUChar3D::IndexType outputSeed12;
		typename TImageUChar3D::IndexType outputSeed21;
		typename TImageUChar3D::IndexType outputSeed22;
		typename TImageUChar3D::IndexType outputSeed23;
		typename TImageUChar3D::IndexType outputSeed24;

		seed[0] = seedlist1[0].toUInt();
		seed[1] = seedlist1[1].toUInt();
		seed[2] = seedlist1[2].toUInt();

		cout << "---- DIVIDE THE LIVER INTO " << m_ui.m_spinParts->value() << " PARTS ----" << endl;

		divideLiverBySeed<TImageUChar3D>(1,
			seed, 
			seed,
			castFilter1->GetOutput(), 
			castFilter2->GetOutput(),
			&outputSk11, 
			&outputSk12,
			&outputLv11,
			&outputLv12,
			&outputSeed11,
			&outputSeed12);

		cout << "outputseed1: " << outputSeed11[0] << ", " << outputSeed11[1] << ", " << outputSeed11[2] << endl;
		cout << "outputseed2: " << outputSeed12[0] << ", " << outputSeed12[1] << ", " << outputSeed12[2] << endl;

		evaluateSeed<TImageUChar3D>(&outputSeed11);
		evaluateSeed<TImageUChar3D>(&outputSeed12);

		if ( m_ui.m_spinParts->value() >= 4)
		{
			// level 2
			divideLiverBySeed<TImageUChar3D>(1,
				outputSeed11, 
				outputSeed11, 
				outputSk11,
				outputLv11,
				&outputSk21,
				&outputSk22, 
				&outputLv21,
				&outputLv22,
				&outputSeed21,
				&outputSeed22);

			divideLiverBySeed<TImageUChar3D>(1,
				outputSeed12, 
				outputSeed12, 
				outputSk12,
				outputLv12,
				&outputSk23,
				&outputSk24, 
				&outputLv23,
				&outputLv24,
				&outputSeed23,
				&outputSeed24);

			cout << "outputseed3: " << outputSeed21[0] << ", " << outputSeed21[1] << ", " << outputSeed21[2] << endl;
			cout << "outputseed4: " << outputSeed22[0] << ", " << outputSeed22[1] << ", " << outputSeed22[2] << endl;
			cout << "outputseed5: " << outputSeed23[0] << ", " << outputSeed23[1] << ", " << outputSeed23[2] << endl;
			cout << "outputseed6: " << outputSeed24[0] << ", " << outputSeed24[1] << ", " << outputSeed24[2] << endl;

			//level 3
			if ( m_ui.m_spinParts->value() >= 8)
			{

				divideLiverBySeed<TImageUChar3D>(1,
					outputSeed21, 
					outputSeed21, 
					outputSk21,
					outputLv21,
					NULL, NULL, NULL, NULL, NULL, NULL);

				divideLiverBySeed<TImageUChar3D>(1,
					outputSeed22, 
					outputSeed22,
					outputSk22,
					outputLv22,
					NULL, NULL, NULL, NULL, NULL, NULL);

				divideLiverBySeed<TImageUChar3D>(1,
					outputSeed23, 
					outputSeed23, 
					outputSk23,
					outputLv23,
					NULL, NULL, NULL, NULL, NULL, NULL);

				divideLiverBySeed<TImageUChar3D>(1,
					outputSeed24, 
					outputSeed24,
					outputSk24,
					outputLv24,
					NULL, NULL, NULL, NULL, NULL, NULL);
			}
		}

		return;

		/*

		TImageType::Pointer m_volumeSkeleton = ((VolumeSource<TImageType>*)m_sourceSkeleton)->getVolume();
		TImageType::Pointer m_volumeLiver = ((VolumeSource<TImageType>*)m_sourceLiver)->getVolume();
		TImageUChar3D::Pointer volumeSkeleton;
		TImageUChar3D::Pointer volumeLiver;


		itk::TimeProbe clock; clock.Start();

		CastImageFilter::Pointer castFilter = CastImageFilter::New();
		castFilter->SetInput(m_volumeLiver);
		castFilter->Update();
		
		ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New();
		maskSO->SetImage(castFilter->GetOutput());
		maskSO->Update();

		TImageUChar3D::RegionType regionType = maskSO->GetAxisAlignedBoundingBoxRegion();
		RegionOfInterestImageFilter::Pointer roiLiverFilter = RegionOfInterestImageFilter::New();
		roiLiverFilter->SetInput(castFilter->GetOutput());
		roiLiverFilter->SetRegionOfInterest(regionType);
		roiLiverFilter->Update();
		
		RegionOfInterestImageFilter::Pointer roiSKFilter = RegionOfInterestImageFilter::New();
		castFilter->SetInput(m_volumeSkeleton);
		castFilter->Update();
		maskSO->SetImage(castFilter->GetOutput());
		maskSO->Update();
		roiSKFilter->SetInput(castFilter->GetOutput());
		roiSKFilter->SetRegionOfInterest(regionType);
		roiSKFilter->Update();

		volumeLiver = roiLiverFilter->GetOutput();		// general pointer
		volumeSkeleton = roiSKFilter->GetOutput();		// general pointer

		TImageUChar3D::IndexType branchIndex;
		TImageType::SpacingType spacing = m_volumeLiver->GetSpacing();

		QStringList seedlist = m_pickingCallBackForVesselSeparation->getPickedPosition().split(" ");
		if ( seedlist.size() != 4)
		{
			cout << "didn't pick any point" << endl;
			return;
		}


		branchIndex[0] = (int(seedlist[0].toDouble()/spacing[0]))-regionType.GetIndex()[0];
		branchIndex[1] = (int(seedlist[1].toDouble()/spacing[1]))-regionType.GetIndex()[1];
		branchIndex[2] = (int(seedlist[2].toDouble()/spacing[2]))-regionType.GetIndex()[2];

		TImageUChar3D::IndexType tmpIndex;

		//TImageType::Pointer m_dupImage = dup1->GetOutput();

		QList<TImageUChar3D::IndexType> branchList;
		QList<TImageUChar3D::Pointer> branchVolumeList;

		for ( int i = branchIndex[0]-1; i <= branchIndex[0]+1; i++)
		{
			for ( int j = branchIndex[1]-1; j <= branchIndex[1]+1; j++ )
			{
				for ( int k = branchIndex[2]-1; k <= branchIndex[2]+1; k++)
				{
					tmpIndex[0] = i;
					tmpIndex[1] = j;
					tmpIndex[2] = k;

					if ( i == branchIndex[0] && j == branchIndex[1] && k == branchIndex[2])
						continue;

					if ( volumeSkeleton->GetPixel(tmpIndex) != 0)
					{
						branchList.push_back(tmpIndex);

						ImageDuplicator::Pointer dup = ImageDuplicator::New();
						dup->SetInputImage(volumeSkeleton); 
						dup->Update();
						dup->GetOutput()->FillBuffer(0);
						dup->GetOutput()->SetPixel(tmpIndex, 1);
						dup->GetOutput()->SetPixel(branchIndex, 1);
						branchVolumeList.push_back(dup->GetOutput());
					}
				}
			}
		}

		for ( int i = branchIndex[0]-1; i <= branchIndex[0]+1; i++)
		{
			for ( int j = branchIndex[1]-1; j <= branchIndex[1]+1; j++ )
			{
				for ( int k = branchIndex[2]-1; k <= branchIndex[2]+1; k++)
				{
					tmpIndex[0] = i;
					tmpIndex[1] = j;
					tmpIndex[2] = k;

					QList<TImageUChar3D::Pointer>::Iterator it = branchVolumeList.begin();
					for (;it != branchVolumeList.end(); it++)
					{
						(*it)->SetPixel(tmpIndex, 1);
					}
				}
			}
		}

		std::cout << "Branch volumes: " << branchVolumeList.size() << std::endl;
		clock.Stop(); printf("(time1: %.2fs)\n", clock.GetMeanTime());
		clock.Start();

		QStack<TImageUChar3D::IndexType> indexStack;
		QList<QList<TImageUChar3D::IndexType>*> indexLists;

		for (int branch = 0; branch < branchList.size(); branch++)
		{
			indexStack.clear();
			indexStack.push(branchList[branch]);
			QList<TImageUChar3D::IndexType>* indexlist = new QList<TImageUChar3D::IndexType>();
			indexLists.push_back(indexlist);

			TImageUChar3D::IndexType nextIndex;
			TImageUChar3D::Pointer dupVolume = branchVolumeList[branch];

			while (1)
			{
				if (indexStack.size() == 0)
					break;

				nextIndex =	indexStack.pop();

				int radius = 1;
				for ( int i = nextIndex[0]-radius; i <= nextIndex[0]+radius; i++)
				{
					for ( int j = nextIndex[1]-radius; j <= nextIndex[1]+radius; j++ )
					{
						for ( int k = nextIndex[2]-radius; k <= nextIndex[2]+radius; k++)
						{
							tmpIndex[0] = i;
							tmpIndex[1] = j;
							tmpIndex[2] = k;

							if ( volumeSkeleton->GetPixel(tmpIndex) !=0 && dupVolume->GetPixel(tmpIndex)==0)
							{
								dupVolume->SetPixel(tmpIndex, 1);

								indexStack.push(tmpIndex);
								indexlist->push_back(tmpIndex);
							}
						}
					}
				}
			}
		}


		QList<QList<TImageType::IndexType>*>::Iterator it = indexLists.begin();
		QList<TImageUChar3D::Pointer>::Iterator vit = branchVolumeList.begin();
		int i = 0;
		for (; it != indexLists.end(); i++)
		{
			if ( (*it)->size() < 100)
			{
				cout << "remove indexlist has " << (*it)->size() << "points" << endl;

				delete *it;
				vit = branchVolumeList.erase(vit);
				it = indexLists.erase(it);
			}
			else
			{
				it++; vit++;
			}
		}

		if ( indexLists.size() != 2)
		{
			cout << "indexlist.size() != 2" << endl;
			return;
		}
		else
		{

			VolumeSource<TImageUChar3D> *subsource1 = 
				new VolumeSource<TImageUChar3D>(QString("Branch#%1").arg(m_sourceSkeleton->rowCount()+1));
			subsource1->setVolume(branchVolumeList[0]);
			m_sourceSkeleton->addChild(subsource1);
			subsource1 = new VolumeSource<TImageUChar3D>(QString("Branch#%1").arg(m_sourceSkeleton->rowCount()+1));
			subsource1->setVolume(branchVolumeList[1]);
			m_sourceSkeleton->addChild(subsource1);
		}

		clock.Stop(); printf("(time 2: %.2fs)\n", clock.GetMeanTime());
		clock.Start();

		SignedDanielssonDistanceMapImageFilter::Pointer dmFilter1 = SignedDanielssonDistanceMapImageFilter::New();
		dmFilter1->SetInput(branchVolumeList[0]);
		dmFilter1->Update();

		SignedDanielssonDistanceMapImageFilter::Pointer dmFilter2 = SignedDanielssonDistanceMapImageFilter::New();
		dmFilter2->SetInput(branchVolumeList[1]);
		dmFilter2->Update();

		clock.Stop(); printf("(time 3: %.2fs)\n", clock.GetMeanTime());

		VolumeSource<TImageUChar3D> *subsource = 
			new VolumeSource<TImageUChar3D>(QString("DistanceMap#%1").arg(m_sourceSkeleton->rowCount()+1));
		subsource->setVolume(dmFilter1->GetOutput());
		m_sourceSkeleton->addChild(subsource);

		VolumeSource<TImageUChar3D> *subsource1 = 
			new VolumeSource<TImageUChar3D>(QString("DistanceMap#%1").arg(m_sourceSkeleton->rowCount()+1));
		subsource1->setVolume(dmFilter2->GetOutput());
		m_sourceSkeleton->addChild(subsource1);

		TImageUChar3D::SizeType size = branchVolumeList[0]->GetBufferedRegion().GetSize();

		cout << "branch volume size [" << size[0] << ", "
									   << size[1] << ", "
									   << size[2] << "]" << endl;


		clock.Start();

		ImageDuplicator::Pointer dup_leftLiver = ImageDuplicator::New();
		dup_leftLiver->SetInputImage(branchVolumeList[0]); 
		dup_leftLiver->Update();
		dup_leftLiver->GetOutput()->FillBuffer(0);

		ImageDuplicator::Pointer dup_rightLiver = ImageDuplicator::New();
		dup_rightLiver->SetInputImage(branchVolumeList[1]); 
		dup_rightLiver->Update();
		dup_rightLiver->GetOutput()->FillBuffer(0);

		TImageUChar3D::Pointer volumeLeftLiver = dup_leftLiver->GetOutput();
		TImageUChar3D::Pointer volumeRightLiver = dup_rightLiver->GetOutput();
		TImageUChar3D::Pointer dmLiver1 = dmFilter1->GetOutput();
		TImageUChar3D::Pointer dmLiver2 = dmFilter2->GetOutput();

		TImageUChar3D::IndexType indexType;
		for ( int i = 0; i < size[0]; i++ )
		{
			for ( int j = 0; j < size[1]; j++)
			{
				for ( int k = 0; k < size[2]; k++)
				{
					indexType[0] = i;
					indexType[1] = j;
					indexType[2] = k;

					if ( volumeLiver->GetPixel(indexType) == 0)
						continue;

					if( dmLiver1->GetPixel(indexType) < dmLiver2->GetPixel(indexType))
					{
						volumeLeftLiver->SetPixel(indexType, 1);
					}
					else
					{
						volumeRightLiver->SetPixel(indexType, 1);
					}
				}
			}
		}

		clock.Stop(); printf("(time: %.2fs)\n", clock.GetMeanTime());

		VolumeSource<TImageUChar3D> *subsourceLeft = 
			new VolumeSource<TImageUChar3D>(QString("LeftLiver#%1").arg(m_sourceSkeleton->rowCount()+1));
		subsourceLeft->setVolume(volumeLeftLiver);
		m_sourceLiver->addChild(subsourceLeft);

		VolumeSource<TImageUChar3D> *subsourceRight = 
			new VolumeSource<TImageUChar3D>(QString("RightLiver#%1").arg(m_sourceSkeleton->rowCount()+1));
		subsourceRight->setVolume(volumeRightLiver);
		m_sourceLiver->addChild(subsourceRight);

		return;*/

		// divide liver

// 
// // 		TImageType::Pointer volume0 = branchVolumeList[0];
// // 		TImageType::Pointer volume1 = branchVolumeList[1];
// 
// 		QList<TImageType::IndexType> indexList0 = *indexLists[0];
// 		QList<TImageType::IndexType> indexList1 = *indexLists[1];
// 
// 
// 		for ( int i = 0; i < size[0]; i++ )
// 		{
// 			std::cout << "[" << i << "]" << std::endl;
// 			for ( int j = 0; j < size[1]; j++)
// 			{
// 				for ( int k = 0; k < size[2]; k++)
// 				{
// 					liverIndex[0] = i;
// 					liverIndex[1] = j;
// 					liverIndex[2] = k;
// 
// 					if ( m_volumeLiver->GetPixel(liverIndex)!= 0 &&
// 						volumeLeftLiver->GetPixel(liverIndex) == 0 &&
// 						volumeRightLiver->GetPixel(liverIndex) == 0)	// construct a cude
// 					{
// 						//TImageType::IndexType tmpIndex;
// 
// 						double minDist0 = 10000000;
// 						double minDist1 = 10000000;
// 
// 						TImageType::IndexType minIndex0;
// 						TImageType::IndexType minIndex1;
// 
// 						double tmp;
// 						for (int l = 0; l < indexList0.size(); l+=13)
// 						{
// 							tmp = (indexList0[l][0] - liverIndex[0])*(indexList0[l][0] - liverIndex[0]) +
// 								(indexList0[l][1] - liverIndex[1])*(indexList0[l][1] - liverIndex[1]) +
// 								(indexList0[l][2] - liverIndex[2])*(indexList0[l][2] - liverIndex[2]);
// 
// 							if ( tmp < minDist0)
// 							{
// 								minDist0 = tmp;
// 								minIndex0[0] = indexList0[l][0]; minIndex0[1] = indexList0[l][1]; minIndex0[2] = indexList0[l][2];
// 							}
// 						}
// 
// 						for (int l = 0; l < indexList1.size(); l+=13)
// 						{
// 							tmp = (indexList1[l][0] - liverIndex[0])*(indexList1[l][0] - liverIndex[0]) +
// 								(indexList1[l][1] - liverIndex[1])*(indexList1[l][1] - liverIndex[1]) +
// 								(indexList1[l][2] - liverIndex[2])*(indexList1[l][2] - liverIndex[2]);
// 
// 							if ( tmp < minDist0)
// 							{
// 								minDist1 = tmp;
// 								minIndex1[0] = indexList1[l][0]; minIndex1[1] = indexList1[l][1]; minIndex1[2] = indexList1[l][2];
// 							}
// 						}
// 
// 						TImageType::IndexType tmpIndex;
// 
// 						for (int l = liverIndex[0]-5; l <= liverIndex[0]+5; l++)
// 						{
// 							for (int m = liverIndex[1]-3; m <= liverIndex[1]+3; m++)
// 							{
// 								for (int n = liverIndex[2]-3; n <= liverIndex[2]+3; n++)
// 								{
// 									tmpIndex[0]=l; tmpIndex[1]=m; tmpIndex[2]=n;
// 
// 									if ( m_volumeLiver->GetPixel(tmpIndex)==0)
// 										continue;
// 
// 									if ( minDist0 < minDist1)
// 									{
// 										volumeLeftLiver->SetPixel(tmpIndex, 1);
// 									}
// 									else
// 									{
// 										volumeRightLiver->SetPixel(tmpIndex, 1);
// 									}
// 								}
// 							}
// 						}
// 					}
// 				}
// 			}
// 		}
// 
// 		VolumeSource<TImageType> *subsourceLeft = 
// 			new VolumeSource<TImageType>(QString("LeftLiver#%1").arg(m_sourceSkeleton->rowCount()+1));
// 		subsourceLeft->setVolume(volumeLeftLiver);
// 		m_sourceLiver->addChild(subsourceLeft);
// 
// 		VolumeSource<TImageType> *subsourceRight = 
// 			new VolumeSource<TImageType>(QString("RightLiver#%1").arg(m_sourceSkeleton->rowCount()+1));
// 		subsourceRight->setVolume(volumeRightLiver);
// 		m_sourceLiver->addChild(subsourceRight);
	}

};



#endif

